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1  Summary 

This  research  program  has  focused  on  enhancing  the  understanding  of  jet  atomization  processes  and  their 
contribution  to  combustion  instabilities  in  liquid  rocket  engines.  During  the  past  three  years,  progress  has 
been  made  in  understanding  the  role  of  the  atomization  process  in  tangential-mode  instabilities  in  the  F- 

1  engine  through  the  use  of  Boundary  Element  Method  (BEM)  codes  developed  with  previous  AFOSR 
sponsorship.  Another  important  research  discovery  resulted  in  the  completion  of  a  zonal  model  capable  of 
addressing  viscous  forces  at  the  periphery  of  the  high-speed  jet.  Through  the  use  of  this  model,  a  jet  swelling 
phenomenon  has  been  fully  characterized;  its  influence  on  the  atomization  process  is  believed  to  be  profound. 

In  addition,  two  new  major  code  developments  have  been  initiated  during  the  subject  contract.  A 
fully  three-dimensional,  nonlinear,  unsteady  BEM  code  has  been  developed  and  validated  for  use  in  future 
atomization  simulations:  This  development  represents  a  major  milestone  in  that  “real  world”  atomization 
processes  are  highly  three  dimensional. 

A  second  new  development,  initiated  in  the  last  few  months  of  the  studies,  involves  a  totally  new  approach 
to  simulating  atomization  problems.  Through  the  use  of  a  fictitious  fluid  “pseudo-density”  which  reflects 
the  local  void  fraction  in  the  two-phase  mixture,  very  complex  atomization  problems  can  be  tackled.  In 
the  model  currently  under  development,  we  project  a  capability  to  track  over  100,000  individual  droplets. 
In  addition,  each  drop  is  subject  to  deformation  from  the  gas-phase  dynamic  pressure  forces,  to  secondary 
atomization,  and  to  possible  collisions  with  other  droplets.  Data  for  these  basic  subprocesses  is  derived  from 
research  results  directly  obtained  by  others  supported  in  AFOSR’s  propulsion  research  program. 

2  Research  Objectives 

The  understanding  of  the  complex  combustion  phenomena  present  in  liquid  rocket  engines  begins  with  the 
fundamental  process  of  fuel  and  oxidizer  jet  atomization.  The  objective  of  this  research  has  been  to  develop 
a  series  of  models,  incorporating  increasingly  complex  physics,  to  assess  the  role  of  atomization  in  the 
combustion  instability  process.  The  models  have  centered  on  the  use  of  Boundary  Element  Methods  (BEMs) 
as  a  means  to  provide  accurate  description  of  these  complex,  nonlinear  processes  under  arbitrary  unsteady 
conditions.  The  models  have  demonstrated  a  capability  to  have  calculations  proceed  beyond  atomization 
events. 

While  the  basic  BEM  techniques  are  inviscid,  recent  development  of  a  zonal  model  using  an  integral 
method  for  boundary  layer  modeling,  permits  a  full  viscous  capability.  This  model,  described  in  Appendix 
A  of  this  report,  is  the  first  primary  atomization  model  to  provide  accurate,  fully  nonlinear  treatment  of 
atomization  processes  under  full-scale  Reynolds  numbers  consistent  with  actual  engine  conditions.  While 
these  BEM  simulations  have  been  useful  in  describing  parent  surfaces  of  modest  complexity,  other  techniques 
are  required  to  resolve  dense  sprays  formed  in  many  rocket  injection  processes.  For  this  reason,  we  have 
embarked  on  the  development  of  a  viscous,  unsteady,  nonlinear  model  capable  of  addressing  flows  in  which 
large  numbers  of  droplets  are  present. 

3  Status  of  the  Research 

Four  major  tasks  have  been  accomplished  in  the  research  project  over  the  past  three  years.  A  viscous 
BEM  simulation  of  a  high-speed  jet,  such  as  those  employed  in  rocket  injectors,  has  been  completed  and 
documented  in  the  open  literature.  This  study  is  summarized  in  Section  3.1.  Simulations  pertaining  to  the 
stability  of  the  F-l  injector  have  been  completed;  these  results  are  summarized  in  Section  3.2.  In  Section 
3.3  we  highlight  some  experimental  results  obtained  in  the  hybrid  rocket  injection  analysis,  while  Section  3.4 
provides  a  insight  into  the  new  3-D  BEM  code  which  has  recently  been  completed.  Finally,  we  conclude  the 
research  status  by  providing  a  description  and  sample  results  from  the  new  dense  spray  model  which  is  still 
under  development. 
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Table  1:  F-l  Injector  Design  Data:  Fuel  Orifice  Test  Data 


Injector  Designation 

Orifice 
Diam.  mm 

Ug,  HZ 

A  P/Pc  % 

Wg/wn 

5U-Flatface 

3.66 

538 

150 

3.0 

5U-Baffled 

4.04 

460 

65 

2.9 

Double  Row  Cluster  (DRC) 

2.79 

454 

400 

1.7 

Prelim.  Flight  Rating  Tests  (PFRT) 

5.79 

- 

0 

5.1* 

Flight  Rating  Test  (FRT) 

7.14 

- 

0 

7.1* 

*  Assumes  ujg  =  460  Hz  IT  chamber  mode. 


3.1  Simulation  of  High-Speed  Viscous  Jets 

An  unsteady,  axisymmetric  BEM-based  model  was  developed.  The  model  included  an  unsteady  integral 
treatment  of  the  thin  boundary  layer  near  the  liquid/wall  and  gas/liquid  boundary.  The  most  remarkable 
aspect  of  these  solutions  is  that  they  all  exhibit  a  “swelling”  in  the  region  very  near  the  orifice  exit  plane. 
Physically,  the  swelling  is  caused  by  the  fluids  transition  from  a  “no-slip”  boundary  condition  inside  the 
orifice,  to  a  free-surface  boundary  condition  outside  the  nozzle.  Figure  1  provides  a  schematic  representation 
of  the  predicted  jet  surface  shape  in  the  region  near  the  orifice  exit.  We  have  found  that  the  amount  of 
swelling  is  dependent  on  the  thickness  of  the  boundary  layer  at  the  orifice  exit  plane  as  illustrated  in  Fig. 
2.  Experimental  verification  of  this  process  is  very  difficult  because  it  represents  a  small  change  (only  a 
few  percent)  in  jet  diameter  in  a  region  which  is  normally  obscured  by  droplets.  However,  we  are  currently 
working  on  indirect  verification  of  the  results  by  correlating  our  predictions  with  jet  cone  angles,  breakup 
lengths,  and  other  measurable  parameters. 

These  results  have  improved  our  understanding  of  the  atomization  process  and  how  the  current  theories 
illustrated  in  Fig.  1  contribute  to  this  process.  In  the  presence  of  a  gas-phase,  the  swelling  of  the  jet  provides 
a  large-amplitude  disturbance  for  initiation  of  aerodynamic  instabilities  which  are  a  critical  element  of  the 
gas  interaction  theory.  For  the  high-speed  jet,  these  instabilities  could  be  amplified  so  fast  that  atomization 
occurs  in  the  region  very  near  the  exit  plane  as  is  observed  in  experiments  under  these  conditions.  Since 
the  amount  of  swelling  is  dependent  on  the  thickness  of  the  boundary  layer  at  the  orifice  exit,  we  can  affect 
the  breakup  characteristics  of  the  jet  by  varying  this  parameter.  For  example,  if  we  wish  the  jet  to  remain 
intact,  we  can  design  an  orifice  passage  to  minimize  the  thickness  of  the  boundary  layer  at  its  exit.  These 
principles  are  used  in  the  design  of  orifices  for  water-jet  cutters  and  fire  hoses.  On  the  other  hand,  if  we 
desire  improved  atomization  (such  as  in  many  combustion  applications),  we  simply  evaluate  designs  which 
promote  thick  boundary  layers  at  the  orifice  exit  plane. 

These  findings  will  help  guide  the  historically  empirical  orifice  design  process  for  combustion  devices.  By 
correlating  performance  of  current  injectors  with  the  swelling  parameter  described  above,  we  will  be  able  to 
tailor  future  designs  to  meet  requirements  without  the  need  of  extensive  (and  expensive)  development  testing. 
Finally,  the  capability  to  predict  atomization  performance  will  lead  to  designs  with  increased  efficiency, 
thereby  not  only  reducing  product  cost,  but  increasing  its  capabilities  at  the  same  time.  Appendix  A 
provides  a  complete  description  of  the  model  and  the  results  of  numerous  simulations. 

3.2  Simulation  of  F-l  Engine  Tangential-Mode  Instability 

In  our  final  report  from  the  previous  contract,  we  provided  some  findings  on  the  stability  of  various  F-l 
injector  designs  used  in  the  engine  development  program.  Here,  2-D  BEM  simulations  were  conducted  in 
which  a  fuel  jet  is  exposed  to  transverse  acoustic  waves  presumed  to  be  present  in  the  chamber.  Both  stable 
and  unstable  injector  designs  from  the  F-l  engine  test  program  have  been  investigated.  As  one  might  expect, 
the  response  of  the  column  grows  dramatically  when  the  acoustic  frequency  (cog)  is  near  that  of  the  column 
natural  frequency  (u>).  Table  1  summarizes  critical  data  for  several  F-l  injector  designs. 

Using  the  model,  a  simulation  was  conducted  of  a  highly-unstable  F-l  injector  configuration  (the  Double 
Row  Cluster,  DRC),  as  well  as  the  final,  stable  configuration  demonstrated  in  Flight  Rating  Tests  (FRT). 
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Orifice's  Wall 


Orifice's  centerline 


Figure  1:  Results  of  computer  simulations  ( which  neglect  the  presence  of  the  gas  phase)  indicate  a  steady, 
non-atomizing  jet  profile.  The  jet  expands,  or  “ swells ”  a  small  amount  near  the  orifice  exit  plane. 


Figure  2:  Swelling  has  been  found  to  be  a  function  of  the  boundary  layer  thickness  at  the  orifice  exit  plane. 
Here,  both  axes  are  plotted  as  a  percentage  of  the  orifice  radius. 
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Figure  3:  Simulation  of  Fuel-Jet  Response  in  the  Unstable  Double-Row  Cluster  (DRC)  and  Stable  Flight 
Rating  Tests  (FRT)  F-l  Engine  Injectors  (times  are  measured  in  seconds  from  start  of  disturbance) 


The  DRC  design  exhibited  a  IT  instability  at  454  Hz  which  led  to  chamber  pressure  oscillations  of  the  order 
of  400%  of  the  mean.  The  most  prominent  difference  between  these  two  injectors  is  the  fuel  orifice  size  (3.57 
mm  radius  on  FRT  vs.  1.4  mm  radius  on  DRC).  The  combination  of  the  acoustic  frequency  in  the  chamber 
and  the  natural  frequency  of  the  DRC  fuel  column  leads  to  conditions  very  near  the  resonant  frequency.  In 
fact,  we  calculate  that  ajg/uj  —  1.7  for  the  DRC,  while  the  FRT  design  has  =  7.0,  a  value  far  from  the 
high  response  region. 

To  assess  the  impact  of  the  fuel  orifice  design  differences  between  the  two  injectors,  we  have  completed  a 
simulation  for  both  designs  at  a  fixed  Weber  number  of  0.1.  Results  of  the  column  shapes  at  various  times 
are  shown  in  Fig.  3  for  both  designs.  This  figure  shows  that  the  unstable  DRC  design  undergoes  violent 
oscillations,  whereas  the  stable  FRT  design  is  relatively  unaffected  by  the  imposed  oscillation.  Clearly,  the 
large  deformations  of  the  DRC  design  will  have  substantial  effect  on  the  jet  impingement  region  -  a  critical 
design  feature  for  this  impinging  element  injector.  For  this  reason,  we  believe  that  the  sensitivity  of 
the  DRC  design  to  transverse  acoustic  energy  may  be  a  major  contributer  to  the  instability 
of  this  injector  design. 

Figure  4  provides  results  for  this  simple  correlation  for  the  five  injector  designs  highlighted  in  Table  1. 
Here,  results  indicate  stability  if  the  frequency  ratio  exceeds  a  value  of  about  3;  beyond  this  value,  the 
column  frequency  is  too  low  to  generate  appreciable  response  from  the  imposed  external  perturbation.  We 
feel  that  these  results  point  to  an  important  new  injector  design  criteria;  the  stability  of  injector  fluid  columns 
should  be  assessed  with  regard  to  coupling  of  the  column  harmonic  frequencies  and  the  acoustic  modes  of 
the  chamber. 

Further  description  of  this  analysis  can  be  found  in  the  Journal  of  Propulsion  and  Power  article  attached 
as  Appendix  B. 
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Frequency  Ratio,  cog  /con 

Figure  4:  Correlation  of  F-l  Engine  Response  with  Frequency  Ratio,  wg/u;n. 

3.3  Hybrid  Rocket  Injection  Experimental  Research 

While  the  main  focus  of  the  studies  in  this  project  is  analytic/numerical,  a  limited  amount  of  experimental 
work  was  accomplished  with  funding  from  this  contract  (as  obtained  through  AFRL  at  Edwards  AFB).  In 
particular,  stable,  reproduceable  ignition  was  verified  in  a  hybrid  rocket  using  hydrogen  peroxide  as  oxidizer 
and  polyethylene  as  fuel.  The  engine  made  use  of  a  unique  Consumable  Catalytic  Bed  (CCB)  ignition  system 
as  highlighted  in  Fig.  5.  The  CCB  permits  rapid,  reproduceable  ignition  events  with  a  minimum  cost  and 
system  complexity;  engines  were  literally  “slam  started”  with  no  special  ignition  valve  opening  procedure. 
A  patent  application  is  in  process  for  the  CCB  design. 

To  investigate  the  functionality  of  the  CCB,  a  transparent  (Plexiglas)  chamber  was  manufactured.  Figure 
6  provides  a  sequence  of  photographs  of  an  actual  test,  showing  the  consumption  of  the  device  and  its  overall 
effect  on  the  flowfield.  Small  perturbations  in  the  pressure  trace  during  some  of  the  tests  are  attributed 
to  ejection  of  small  portions  of  CCB  during  the  end  stages  of  consumption.  Overall,  this  injection/ignition 
concept  is  very  attractive  for  hybrid  engines  not  requiring  a  restart  capability.  Further  details  regarding  this 
project  are  summarized  in  the  manuscript  in  Appendix  C. 

3.4  Development  of  a  3-D  BEM  Code 

A  substantial  new  development  resulting  from  this  research  involves  the  creation  of  a  full  three  dimensional 
computational  tool.  The  model  employs  linear  elements  and  has  undergone  extensive  validation  on  both 
droplet  and  liquid  jet  problems.  Sample  liquid  jet  grids  are  shown  in  Fig.  9  for  various  node  densities. 
The  code  was  run  on  a  classic  “infinite  jet”  simulation  in  order  to  assess  performance  against  analytic 
results.  Figure  10  summarizes  the  results  of  the  calculations;  the  growth  rate  (surface  deformation  rate) 
is  compared  against  the  analytic  result  for  a  variety  of  wavenumbers,  k .  At  low  wavenumbers,  excellent 
agreement  is  obtained.  However,  the  results  at  high  wavenumbers  are  not  as  accurate  due  to  the  fact  that 
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Figure  5:  Cross  Sectional  View  of  the  Combustion  Chamber 


the  computational  meshes  employed  are  not  sufficient  to  resolve  the  surface  shape  to  the  desired  accuracy. 
These  computations  turned  out  to  be  very  expensive;  a  single  simulation  here  could  take  as  long  as  one  week 
on  a  state-of-the-art  computer.  More  work  needs  to  be  performed  to  speed  up  the  code  through  the  use  of 
parallelization  and  advanced  matrix  inversion  schemes. 

Another  issue  which  warrants  discussion  is  the  surface  treatment.  Originally,  we  had  arranged  with 
Professor  Chandrajit  Bajaj,  a  leader  in  scientific  visualization  with  a  large  research  group  in  Computer 
Science,  to  provide  surface  remeshing  capabilities  for  our  simulations.  Unfortunately,  Dr.  Bajaj  left  Purdue 
for  a  position  at  the  University  of  Texas  in  early  1998.  Prior  to  his  departure,  his  students  were  having 
great  difficulty  adapting  their  current  remeshing  schemes  to  our  problems  since  we  typically  demand  much 
higher  accuracy  than  that  of  a  simple  3-D  visualization  of  some  object.  For  this  reason,  the  present  model 
does  not  contain  a  remeshing  capability;  a  factor  which  negates  the  use  of  the  code  for  highly  nonlinear 
calculations.  In  these  cases,  mesh  distortion  due  to  surface  motion  makes  regridding  a  necessity.  This  area 
is  one  of  current  research;  we  will  continue  to  explore  mechanisms  to  move  the  surface  in  a  more  accurate 
and  effective  fashion. 

Appendix  D  provides  additional  documentation  regarding  the  3-D  model  development  and  validation. 

3.5  Development  of  a  Dense  Spray  Model 

The  last  task  initiated  in  the  three-year  research  program  involved  the  early  stages  of  development  of  a 
dense  spray  model.  Here,  our  motivation  is  to  provide  a  tool  for  analysis  of  very  complicated  atomization 
problems  (such  as  those  in  liquid  rocket  engines)  in  which  105  —  106  droplets  may  be  present.  A  unique 
“quasi-axisymmetric”  methodology  has  been  developed  in  which  we  analyze  the  spray  produced  along  a  single 
plane  intercepting  the  axis  of  the  jet/atomizer.  Droplets  are  launched  from  the  inner  core  of  the  jet  and 
are  permitted  to  collide  and  undergo  secondary  atomization.  Droplet  collision  and  secondary  atomization 
modules  draw  from  previous  AFOSR  research  on  these  fundamental  processes.  For  example,  we  utilize 
results  from  Professor  Faeth’s  drop  atomization  work  at  Michigan  and  from  Professor  Law's  drop  collision 
work  at  Princeton  -  both  projects  were  supported  by  AFOSR  funding. 

These  droplet  collision/atomization  modules  permit  the  tracking  of  roughly  100,000  droplets  in  a  mas- 


7 


Figure  6:  Still  Photographs  of  Motor  Y1  showing  CCB  consumption.  Here,  the  CCB  is  located  on  the  right 
hand  side  of  the  image  and  flow  is  from  right  to  left.  From  top  to  bottom,  images  are  for  t  =  3.3,  5.45,  5.5, 
and  11.5  sec 
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Figure  7:  Sample  Grids  Utilized  in  3-D  BEM  Liquid  Jet  Simulations 
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Figure  8:  Validation  of  3-D  BEM  on  Linear  Growth  Rates  (u;)  of  a  Liquid  Jet 


sively  parallel  calculation.  The  coupling  of  the  flowfield  with  the  gas-phase  behavior  is  accomplished  by 
a  parallel  unsteady  Navier-Stokes  calculation  which  presumes  the  two-phase  flow  can  be  represented  by  a 
pseudo-fluid  which  varies  in  density  between  the  liquid  and  gas  extremes.  We  have  had  substantial  success 
using  a  pseudo-fluid  implementation  to  simulate  complex  cavitating  flows  in  which  large  numbers  of  bubbles 
are  present.  For  this  reason,  the  methodology  is  also  attractive  for  sprays  in  which  large  numbers  of  droplets 
are  present. 

The  model,  which  is  the  subject  of  a  follow-on  grant  by  AFOSR  and  Dr.  Mitat  Birkan,  is  in  the  final 
stages  of  validation.  In  the  coming  year,  we  should  have  the  capability  to  perform  complete  simulations  of  a 
dense  spray  and  actually  generate  spray  statistics  (such  as  mean  diameter,  mass  fluxes,  and  velocity  fields) 
for  the  entire  spray.  This  development  will  provide  a  tool  useful  for  real  sprays  of  interest  to  both  the  liquid 
rocket  and  atomization  communities.  A  more  detailed  description  of  the  model  is  contained  in  Appendix  E 
of  this  report. 

4  Professional  Activities 

The  efforts  outlined  in  the  previous  section  of  this  report  were  made  possible  by  two  grants  from  AFOSR. 
A  single  student,  Mr.  Chienchi  Chao,  was  supported  under  the  base  grant.  In  addition,  some  funds  were 
obtained  from  AFRL  at  Edwards  AFB;  these  funds  appeared  as  augmentations  to  the  initial  award  for  the 
present  contract.  These  monies  were  used  to  support  hybrid  rocket  injection  studies  and  the  Ph.D.  student 
Dr.  Eric  Wernimont. 

Ph.D.  Dissertations 


Hilbing,  J.  M.,  “Modeling  Liquid  Jet  Atomization  Processes,”  August,  1996. 
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Wernimont,  E.,  “Experimental  Study  of  Combustion  in  Hydrogen  Peroxide  Hybrid  Rockets,”  August, 
1997. 

Chao,  C.,  “Three-Dimensional  Nonlinear  Jet  Atomization  Model,”  December,  1998. 

M.S.  Theses 


Rump,  K.,  “Modeling  the  Effect  of  Unsteady  Chamber  Conditions  on  the  Atomization  Process,”  Decem¬ 
ber,  1996. 

Murray,  I.,  “Modeling  Acoustically  Induced  Oscillations  of  Droplets,”  August,  1996.  (A  paper  involving 
this  work  was  awarded  First  Place  in  the  AIAA  Graduate  Student  Conference  held  in  April,  1996  at  Purdue 
University.) 


A  list  of  journal  publications  (and  submissions)  associated  with  these  efforts  are  provided  in  the  following 
list.  Highlighted  items  (*)  have  been  attached  in  the  Appendices  of  this  report. 

Refereed  Journal  Publications  and  Submissions 


1.  Heister,  S.  D.,  Rutz,  M.,  and  Hilbing,  J.,  “Effect  of  Acoustic  Perturbations  on  Liquid  Jet  Atom¬ 
ization”,  Journal  of  Propulsion  and  Power ,  V13,  No.  1,  pp.  82-88,  1997. 

2.  Heister,  S.  D.,  “Boundary  Element  Methods  for  Two-Fluid  Free  Surface  Flows”,  V19,  No.  4, 
Engineering  Analysis  with  Boundary  Elements ,  pp.  309-317,  1997. 

3.  *Hilbing,  J.  H.,  and  Heister,  S.D.,  “Nonlinear  Simulation  of  a  High-Speed,  Viscous,  Liquid  Jet”, 
Atomization  and  Sprays ,  V  8,  pp.  155-178,  1997. 

4.  Rump,  K.  M.,  and  Heister,  S.  D.,  “Modeling  the  Effect  of  Unsteady  Chamber  Conditions  on 
Atomization  Processes”,  Journal  of  Propulsion  and  Power ,  V  14,  pp.  576-578,  1998. 

5.  Murray,  I.  F.,  and  Heister,  S.  D.,  “On  a  Droplet’s  Response  to  Acoustic  Excitation”,  International 
Journal  of  Multiphase  Flow ,  To  Appear,  1999. 

6.  Wernimont,  E.J.,  and  Heister,  S.D.,  “A  Reconstruction  Technique  for  Reducing  Hybrid  Rocket 
Combustion  Test  Data”,  Journal  of  Propulsion  and  Power ,  To  Appear,  1999. 

7.  *Chao,  C.  and  Heister,  S.  D.,  “Contributions  of  Atomization  to  F-l  Engine  Combustion  Instabil¬ 
ities”,  Journal  of  Propulsion  and  Power ,  In  Review,  1997. 

8.  *Wernimont,  E.J.,  and  Heister,  S.D.,  “Combustion  Experiments  in  a  Hydrogen  Peroxide  Polyethy¬ 
lene  Hybrid  Rocket  with  Catalytic  Ignition”,  In  Review,  Journal  of  Propulsion  and  Power,  1998. 

A  list  of  the  conference  papers  presented  in  association  with  work  under  these  grants  is  provided  in  the 
list  below.  The  starred  item  is  included  in  Appendix  E  of  this  report. 

Conference  Papers  and  Presentations 


1.  Hilbing,  J.H.,  Heister,,  S.D.,  and  Rump,  K.,  “Recent  Advances  in  Nonlinear  Modeling  of  Atom¬ 
ization  Processes”,  ILASS-96  Conference  Proceedings,  1996. 

2.  Murray,  I.  F.,  and  Heister,,  S.D.,  “Modeling  Acoustically-Induced  Oscillations  of  Droplets”, 
ILASS-96  Conference  Proceedings,  1996. 

3.  Wernimont,  E.  J.,  and  Heister,  S.  D.,  “Progress  in  Hydrogen  Peroxide  Oxidized  Hybrid  Rocket 
Experiments”,  AIAA  96-2696,  32nd  AIAA  Joint  Propulsion  Conference,  1996. 
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4.  Rump,  K.  M.,  and  Heister,  S.  D.,  “Modeling  the  Effect  of  Unsteady  Chamber  Conditions  on 
Atomization  Processes”,  AIAA  97-3298,  33rd  AIAA  Joint  Propulsion  Conference,  1997. 

5.  Heister,  S.  D.,  “Modeling  Atomization  Processes”,  52nd  Annual  Industrial  Waste  Conference 
Proceedings,  1997. 

6.  *Heister,  S.  D.,  “Modeling  Primary  Atomization  Processes”,  AIAA  98-3837,  34th  AIAA  Joint 
Propulsion  Conference,  1998. 

7.  Heister,  S.  D.,  Wernimont,  E.  J.,  and  Rusek,  J.  J.,  “High  Test  Peroxide  Hybrid  Rocket  Re¬ 
search”,  Hydrogen  Peroxide  Propulsion  Workshop,  Surrey  England,  July,  1998.  (see  also:  http  : 
/  /  www.ee. surrey. ac.uk/ EE /CSER/U  OS  AT  /  conf  /  sheister.htm) 

4.1  Technology  Transfer/Coupling  Activities 

Numerous  technology  transfers  have  occurred  during  the  period  associated  with  these  grants.  These  items  are 
summarized  the  the  table  provided  above.  Models  currently  under  development  should  be  of  great  interest 
to  the  liquid  and  hybrid  rocket  engine  community. 

TECHNOLOGY  TRANSFER 


Performer 

Customer 

Result 

Application 

1 

S.  D.  Heister 

Purdue  University 
Ph:  (317)  494-5126 

KB  Sciences 

R.  Humble 

Ph:  (719)  333-6554 

Provided  Injector  Designs 

KEW  Engines  using 
Non-Toxic  Propellants 

2 

S.  D.  Heister 

Purdue  University 
Ph:  (317)  494-5126 

Allied  Signal 

B.  Mar  riot 

Ph:  (219)  254-5235 

Ramjet  Inlet  & 

Combustor  Analysis 

Talos  Missile 

Target  Drone 

3 

S.D.  Heister 

Purdue  University 
Ph:  (317)  494-5126 

TRW  Space  Group 

J.  Calvin ac 

Ph:  (310)  812-5314 

Pintle  Injector  Design 

Shuttle  Non-Toxic 
LOX/Ethanol  Thruster 

4 

S.D.  Heister 

Purdue  University 
Ph:  (317)  494-5126 

TRW  Space  Group 
R.  Sackheim 

Ph:  (310)  813-9304 

Pintle  Injectors 

Combustion  Instability 
Studies 

5 

S.  D.  Heister 

Purdue  University 
Ph:  (317)  494-5126 

KB  Sciences 

R.  Humble 

Ph:  (719)  333-6554 

Designed  Injector 

LLNL  GOX/GH2 

Water  Rocket  Concept 
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5  Appendix  A  -  Viscous  Jet  Simulations 


Hilbing,  J.  H.,  and  Heister,  S.D.,  “Nonlinear  Simulation  of  a  High-Speed,  Viscous,  Liquid 
Jet”,  Atomization  and  Sprays,  V  8,  pp.  155-178,  1997. 
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NONLINEAR  SIMULATION  OF  A  HIGH-SPEED,  VISCOUS  LIQUID  JET 
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*  Staff  Engineer,  TRW  Space  &  Electronics  Group,  Redondo  Beach,  CA  90278 

t  Associate  Professor,  School  of  Aeronautics  and  Astronautics,  Purdue  University,  W.  Lafayette,  IN  47907 


Abstract 

A  model  has  been  developed  to  simulate  the  nonlinear,  unsteady  evolution  of  a  high-speed  viscous  liquid 
jet  issuing  from  a  circular  orifice.  The  model  is  based  on  a  zonal  approach  in  which  an  integral  method 
is  utilized  for  a  thin  viscous  region  at  the  jet  periphery,  while  a  boundary  element  method  is  used  for  the 
inviscid  “core”  flow.  Results  indicate  that  steady-state  solutions  are  possible  neglecting  the  presence  of  the 
gas.  Under  these  conditions,  the  jet  “swells”  in  diameter  and  the  boundary  layer  thins  to  a  shear  layer 
over  a  length  of  about  half  an  orifice  radius.  Because  boundary  layer  relaxation  is  occurring  during  these 
simulations  in  which  steady-state  solutions  appear,  atomization  mechanisms  relying  on  this  process  cannot 
explain  the  observed  behavior.  The  swelling  phenomenon  has  the  potential  to  explain  several  fundamental 
experimental  atomization  observations  regarding  turbulence  and  orifice  design. 


Introduction 

In  many  applications  (combustors,  paint,  cosmetic,  and  agricultural  sprays  to  name  a  few),  it  is  highly 
desirable  to  inject  a  liquid  in  such  a  fashion  that  rapid  atomization  takes  place.  In  spite  of  the  high  interest 
in  spray-producing  injectors,  the  design  of  these  devices  is  largely  empirical  due  to  a  lack  of  knowledge  of 
the  basic  atomization  mechanisms.  Understanding  of  these  mechanisms  is  incomplete  because  of  difficulties 
in  making  experimental  observations  very  near  the  nozzle  and  in  modeling  the  often  turbulent,  two-phase 
flow  in  this  region. 

At  the  present  time,  several  mechanisms  have  been  suggested  to  contribute  (or  explain)  the  atomization 
process.  Aerodynamic  interactions  are  of  obvious  importance  in  generating  Kelvin-Helmholtz  type  instabil¬ 
ities  and  in  causing  secondary  atomization  of  larger  drops  shed  from  the  periphery  of  the  jet.  Theory  for 
this  mechanism  is  due  to  G.I.  Taylor1,  with  more  recent  extensions  and  modifications2.  Unfortunately,  this 
mechanism  cannot  explain  all  atomization  behavior.  Several  researchers3’4  have  shown  that  the  design  of  the 
orifice  passage  definitely  contributes  to  atomization  behavior  in  many  instances.  Additionally,  cavitation5,6 
in  the  orifice  also  can  play  an  important  role  in  some  cases. 

Many  researchers  have  proposed  that  the  state  of  the  boundary  layer  exiting  the  orifice  is  a  critical 
atomization  mechanism.  For  example,  some  have  suggested  that  boundary  layer  turbulence  plays  a  critical 
role  in  the  atomization  process.  While  turbulence  may  play  a  role  in  destabilizing  low-speed  jets7,8,  at 
the  higher  speeds  consistent  with  the  atomization  regime,  results  indicate3,9  that  nozzles  with  turbulent 
boundary  layers  provided  poorer  atomization  (as  evidenced  by  a  smaller  cone  angle)  than  laminar  jets.  In 
addition,  theoretical  work10,11  suggests  that  turbulent  velocity  profiles  are  actually  more  stable  than  the 
more  fully-developed  profiles  associated  with  laminar  flows. 

Other  authors12,13  have  suggested  that  boundary  layer  velocity  profile  rearrangement  (due  to  the  change 
from  a  no-slip  to  a  free  surface  condition  at  the  orifice  exit)  can  explain  atomization  behavior.  While  there 
is  some  doubt  that  this  mechanism  can  fully  explain  atomization3,  this  factor  has  not  been  eliminated  as  a 
potential  contributor  at  this  point  in  time. 

Recent  development  of  fully-nonlinear  models  based  on  boundary  element  methods14-16  have  enhanced 
abilities  to  simulate  complex  atomization  processes  under  the  assumption  of  inviscid  flow.  Most  recently, 
models  have  been  created16-18  to  include  the  effects  of  orifice  geometry  on  the  atomization  process.  The 
present  work  provides  an  extension  to  these  previous  models  via  the  incorporation  of  a  viscous  region  at  the 
periphery  of  the  jet  (as  well  as  inside  the  orifice  passage).  The  model  is  applied  over  a  range  of  Reynolds 
and  Weber  numbers  to  shed  light  on  the  role  of  boundary  layer  development  on  the  atomization  process. 
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Model  Development 

At  the  present  time,  there  are  very  few  analytic/numerical  tools  available  to  provide  high- resolution,  nonlinear 
simulations  of  liquid  jets  in  the  atomization  regime.  This  regime  is  typically  characterized  by  thin  boundary 
layers  (typically  less  than  10%  of  the  orifice  radius)  and  high  Reynolds  numbers  (of  the  order  of  104-105). 
While  the  fully- viscous  treatments  of  Unverdi  and  Tryggvason19  and  Osher  and  Sethian20  have  demonstrated 
abilities  to  provide  high- resolution  simulations  of  capillary  flows,  these  methods  are  currently  limited  (due 
to  computational  resources)  to  Reynolds  numbers  of  the  order  of  102. 

Boundary  element  methods  (BEMs)  are  quite  attractive  for  inviscid  problems,  however  the  nonlinearity 
of  the  governing  equations  for  viscous  flow  present  substantial  complications  for  this  modeling  technique. 
Recent  development  of  the  Dual  Reciprocity  Method21  has  extended  BEM  capabilities  to  viscous  flows,  but 
this  technique  also  has  trouble  in  providing  adequate  resolution  for  high  Reynolds  number  problems17. 

For  this  reason,  it  appears  that  a  zonal  approach  which  provides  a  separate  treatment  for  the  viscous 
region  would  be  attractive.  Such  techniques  have  been  supplied  successfully  by  Wu22,23  by  solving  a  set 
of  boundary  layer  equations  (by  standard  CFD  methods)  in  viscous  regions,  and  using  BEM  at  the  edge 
of  the  boundary  layer.  In  free  surface  problems  with  moving  grids,  application  of  a  full  CFD  grid  to  the 
boundary  layer  would  be  a  substantial  expense  in  today’s  computing  environment.  For  this  reason,  we  have 
chosen  an  integral  method  to  model  the  behavior  of  the  viscous  zone.  Hailey,  et.  al.24  had  successfully 
adopted  this  technique  to  solve  viscous  external  flows  over  solid  boundaries.  Using  this  approach  in  an 
axisymmetric  unsteady  calculation,  the  entire  jet  can  be  modeled  with  two  lines  of  nodes:  one  at  the 
gas/liquid  or  liquid/ wall  interface,  and  one  at  the  edge  of  the  liquid  boundary  layer. 

We  presume  that  all  variables  have  been  nondimensionalized  by  choosing  the  liquid  density,  the  orifice 
radius,  and  the  ideal  (Bernoulli)  orifice  exit  velocity  as  dimensions.  When  compared  with  the  well-known 
Karman-Pohlhausen  technique,  the  present  analysis  is  complicated  due  to  unsteadiness  in  the  boundary 
layer,  and  due  to  the  fact  that  the  boundary  layer  is  axisymmetric  (rather  than  two  dimensional).  This 
latter  point  is  important,  since  the  boundary  layer  will  naturally  tend  to  be  reduced  in  thickness  under 
outward  radial  motion  due  to  conservation  of  mass  considerations.  The  governing  equations  for  the  viscous 
region  can  be  developed  by  assuming  a  polynomial  approximation  for  the  velocity  profile,  solving  for  the 
coefficients  of  the  polynomial  based  on  edge  conditions,  and  finally  substituting  the  assumed  profile  into  the 
integral  form  of  the  momentum  equation  and  completing  the  integrations.  The  following  subsections  outline 
this  procedure  for  two  cases:  a  boundary  layer  on  a  solid  wall,  and  a  boundary  layer  on  a  free  surface. 


Solid  Wall  Boundary  Layer 

The  integral  analysis  derived  herein  provides  an  extension  to  the  results  derived  by  Howarth25  for  the  case 
of  a  unsteady,  axisymmetric  flow  over  a  solid  boundary.  Here,  we  have  unique  treatments  due  to  the  fact 
that  we  have  an  internal  flow  within  the  orifice  passage,  and  because  of  the  free-surface  boundary  condition 
applied  outside  the  orifice.  The  boundary  layer  profile  is  assumed  to  be  approximated  by  a  fourth-order 
polynomial: 

—  =  aQ  +  ai  T]  +  a2  if)2  +  a3  rj3  +  a4  rf  (1) 

tig 

where  u  and  ue  are  velocities  parallel  to  the  wall  within  and  at  the  edge  of  the  boundary  layer,  respectively. 
The  parameter  rj  =  y/8  varies  from  0  at  the  wall  to  1  at  the  edge  of  the  inviscid  core.  The  five  coefficients 
(a0  —  04)  are  determined  from  the  following  boundary  conditions26: 

“<°>=teU=(0U=° 

u(S)  =  ue  ;  (fF)y=0  =  Re(H) 

The  fifth  condition  in  Equation  2  comes  from  examination  of  the  ^-momentum  equation  at  the  wall,  using 
u( 0)  =  u(0)  =  0  and  d2u/dx 2  d2u/dy2 .  Solving  for  the  coefficients  in  Equation  1  using  the  given  boundary 
conditions  yields: 

7-  =  (2»?-  2t?3  +  »?4)  +  j  (1  -  rjfv  (3) 
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where 


A  =  _£S?(|f)  (4) 

ue  \ax ) 

is  a  pressure  gradient  parameter,  and  is  derived  from  examination  of  the  unsteady  momentum  equation  at 
the  outer  edge  of  the  boundary  layer. 

The  governing  integral  momentum  equation  can  be  derived  by  considering  an  incremental  width  of  the 
boundary  layer  at  some  radial  location  ro  of  the  orifice  wall,  as  shown  in  Figure  1.  An  x-y  axis  system  is 
aligned  with  the  orifice  wall  such  that  the  y- axis  extends  into  the  fluid  and  makes  an  angle  of  C  with  a  line 
in  the  radial  direction.  Note  that  £  =  0  for  a  simple,  non-tapered  orifice  passage.  A  momentum  balance  in 
this  annular  section  of  the  boundary  layer  follows  the  derivation  given  by  Howarth25,  and  for  this  unsteady, 
internal  flow  can  be  written: 


where 


y  =  r0  -  yc osC 


6*  is  the  displacement  thickness,  and  9  the  momentum  thickness,  defined  as: 


(6) 


^coe<)  0  -  £) dy 
s=[  0  -  tos<)  0  -  z)  (£)ds 


respectively.  Substituting  the  assumed  profile  from  Equation  3  into  Equations  7  and  8  results  in  integrals 
which  can  be  evaluated  analytically.  Performing  these  integrations  gives: 


for  the  displacement  thickness,  and: 


9  =  6 


37  A 
315  +  945 


A2  S  J_ 5_ 

9072  r0  C°S  ^  V 126 


_A _ tf_\ 

945  30240  / 


(9) 


(10) 


for  the  momentum  thickness. 

Equations  9  and  10  can  be  substituted  into  Equation  5  to  give  a  single  equation  involving  8,82,  (08/ Ox), 
and  8(05/ dx).  The  remaining  quantities  in  the  equation  (A,A2,C,  etc.)  are  known  along  the  boundary 
layer  interface  since  we  presume  that  dp/dx  =  -ue  (0ue/0x)  in  accordance  with  standard  boundary  layer 
assumptions.  The  momentum  equation  can  be  solved  by  utilizing  a  first-order  backwards  difference  to 
approximate  the  derivative  of  the  boundary  layer  thickness  at  node  i: 


(11) 


Assuming  that  the  8  and  S2  quantities  in  the  equation  are  evaluated  at  node  i,  the  equation  reduces  to  a 
cubic  equation  for  the  boundary  layer  thickness: 


c3  Sf  +  c2  Sf  +  ci  Si  +  c0  =  0  (12) 

which  can  be  solved  for  rt;.  Since  <J,_i  is  assumed  known  and  is  included  in  the  c,-  coefficients  in  Equation  12, 
this  equation  can  be  solved  by  marching  along  the  boundary  layer  using  an  initial  value  of  the  boundary 
layer  thickness  at  the  first  node.  This  differencing  and  solution  procedure  is  consistent  with  the  parabolic 
nature  of  the  boundary  layer  equations.  The  c,-  values  are  quite  complicated  functions  whose  definitions  are 
included  in  the  Appendix. 
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Free  Surface  Boundary  Layer 

Equation  1  is  again  assumed  to  approximate  the  boundary  layer  profile  along  the  free  surface  of  the  jet. 
Here,  the  boundary  conditions  are  chosen  as: 


u(0)  =  us  ;  u(S)  -  ue 


=  0 


so  that  their  substitution  yields  the  boundary  layer  profile: 


u 

ue 


(i_^)  (6,2_87?3  +  37?4)  +  ^ 


(13) 


(14) 


where  us  is  the  velocity  at  the  surface  of  the  jet. 

Figure  1  shows  the  nomenclature  for  the  boundary  layer  along  the  free  surface  of  the  jet.  As  in  the  case 
of  the  solid  wall,  the  integral  momentum  equation  is  derived  by  applying  conservation  of  momentum  to  the 
annular  region  containing  the  boundary  layer.  This  equation  can  be  written  in  general  form  as: 


iff  •  j  u  dm  +  J  u  dm  =  j  aijiii  dA^ 


(15) 


where  V  and  S  denote  the  volume  and  surface  area  of  the  annular  section,  respectively,  ex  denotes  the 
tangential  unit  vector,  dm  =  p  dVf  and  dm  =  p(u  —  us)h  dA.  Here,  Cij  is  the  stress  tensor,  defined  as: 


-pSij  +  p 


(dui  duj  \ 

dxj  dxi  J 


(16) 


with  Sij ,  the  Kronecker  delta.  Carrying  out  the  dot  product,  the  first  term  can  be  written: 

‘•■{ft  [X 5  H } = s  [e-  ■  X s  *■]-£•  [X  * dm 


which  can  be  further  simplified  to: 


[X  *  H } = s  IX  “  H  ■  (^).  X  ”  *" 


(17) 


(18) 


Note  that  dex/dt  is  nonzero  during  free  surface  motion,  but  only  has  a  component  in  the  direction  normal 
to  the  surface,  denoted  with  a  subscript  y  in  the  equation.  In  evaluating  the  right  hand  side  of  Equation  15, 
second  derivatives  of  the  ^-velocity  are  neglected  in  accordance  with  standard  boundary  layer  theory26. 
However,  the  ^-derivatives  of  the  y-velocity  are  retained.  Using  these  assumptions,  and  Equation  18,  the 
momentum  equation  in  integral  form  becomes18’25: 


where  y  is  defined  as: 

y  =  r0-ycosC  (20) 

Here,  evaluation  of  volume  integrations  leads  to  terms  involving  conditions  at  the  boundary  shared  with 
the  inviscid  (potential)  core.  For  example,  the  term  involving  d<f>/dy  represents  the  velocity  normal  to  the 
inner  edge  of  the  boundary  layer  shown  in  Fig.  1.  The  quantity  <j>  represents  the  velocity  potential  in  the 
inviscid  core  region.  If  the  assumed  boundary  layer  profile  in  Equation  14  is  substituted  into  Equation  19, 
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the  integrals  can  be  evaluated,  and  the  result  is  an  equation  involving  5,  (85/dx),  8(88 /dx),  ( 88/dt )  and 
6(85 /dt).  The  remaining  quantities  are  known  along  the  entire  surface.  In  the  same  manner  used  in  the 
previous  section,  backwards  differences  can  be  used  to  approximate  (85/dx)  and  (88/dt): 

(88\n  S?-S?_ i  f  88 \n  _  8?  -  6?-1 

\8x)i  “  Ax  \8t)i  “  At 

Substitution  of  these  approximations  into  Equation  19  results  in  an  equation  of  the  form: 

bo  +  b\  8/  +  62  (8" )2  =  0  (21) 

where  the  subscript  i  denotes  the  node,  and  the  superscript  n  ,  the  time  level.  Values  of  the  boundary 
layer  thickness  at  previous  nodes  and  previous  time  levels  are  present  in  the  known  60  and  61  coefficients  of 
Equation  21,  so  a  time  and  space  marching  method  can  be  used  to  solve  for  the  boundary  layer  thickness 
along  the  free  surface.  In  the  current  method,  8  is  calculated  starting  at  the  orifice  exit,  and  the  boundary 
layer  in  the  orifice  must  be  calculated  before  that  for  the  free  surface.  The  Appendix  contains  the  rather 
involved  definitions  of  the  6,-  coefficients. 


Surface  Velocities  and  Flow  Kinematics 

Solution  of  Laplace’s  equation  in  the  inviscid  core  results  in  components  of  the  velocity  vector  (8<j>/ds  and  q) 
at  the  inner  edge  of  the  boundary  layer16.  Two  additional  conditions  are  required  to  determine  the  tangential 
and  normal  velocities  at  the  free  surface. 

An  examination  of  the  x-momentum  equation  on  the  free  surface  yields  an  equation  for  the  tangential 
velocity  on  the  surface,  us.  In  Lagrangian  form,  the  dot  product  of  the  vector  momentum  equation  and 
tangential  unit  vector  can  be  written: 


.  (Du  „  1  _j,\ 


(22) 


Carrying  out  the  dot  product  produces  another  term  involving  the  time  derivative  of  the  tangential  unit 

vector:  _  ^  „ 

Dus  _  Dex  dp  t  1  2 


-  t  x  rr2  ZOO\ 

Dt  ~U"Dt~~8i  +  ^  Us  (  3) 

Since  only  the  normal  component  of  this  change  in  the  tangential  unit  vector  is  non-zero,  Equation  23  reduces 

Dus  ( Dex  \  dp  J_  /  d^Us  .  d2us ' 


Dt 


Qs 


Dt 


+ 


dx  Re  \  dx2  dy 2 


(24) 


Again,  ( 82us/8x 2)  is  neglected  in  the  boundary  layer,  and  substitution  of  the  boundary  layer  profile  permits 
evaluation  of  the  remaining  second  derivative: 

i2<y2  (ue  -  us ) 


Dus  _  / Dix\  _  8p_ 

~Dt~qs  \  Dt  )  ~  8x 


+ 


Re 


(25) 


Equation  25  is  integrated  in  time  to  solve  for  the  change  in  tangential  velocity  along  the  free  surface. 

The  normal  velocity  at  the  surface  can  be  derived  by  examination  of  the  integral  continuity  equation  in 
the  annular  section  shown  in  Figure  1.  This  equation  can  be  solved  for  q$  to  yield: 

r5 


Qs  =  -  1 - cos  C 

T  o 


dl 

dy 


_  ]_d_  r 
ro  dx  J0 


yu  dy 


(26) 


where  d</)/dy  is  defined  as  positive  in  the  direction  into  the  fluid.  Substituting  the  boundary  layer  profile 
into  this  equation,  integrating  to  the  boundary  layer  thickness,  and  differentiating  the  last  term  with  respect 
to  a?,  gives: 

A  s  \  d<j> 

„=_h__cos(j_ 


+ 


[( 


5 

5_ 

10 


2cosC  d8 
tq  dx 


£sinC  d( 
ro  dx 


(4  ue  +  us)  + 


6  cos£ 


4 

dx  dx 


(27) 
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For  small  S,  dilitation  and  rotation  of  the  boundary  layer  due  to  changes  in  r0,  ue,  and  u,  is  minor  and  the 
equation  reduces  to  the  straightforward  result  qs  ~  —d<j>/dy.  Inside  the  orifice  passage,  the  solid  surface 
dictates  that  qs  =  0  in  this  region. 

For  the  node  at  the  tip  of  the  jet  (on  the  jet  centerline),  Eq.  27  becomes: 


(?*)r0=0 


dy 


'  due  dus 


d£ 

dx 


+ 


du^ 

dx 


) 


where  we  have  used  the  basic  geometric  result: 


lim 

Tq  — hO 


cosC 

ro 


dC 

dx 


(28) 


(29) 


in  deriving  this  relation. 

The  surface  is  tracked  along  its  instantaneous  velocity  vector  in  a  manner  similar  to  the  method  used  in 
Hilbing,  et.  al.16  with  axial  and  radial  components  of  the  surface  velocity  given  by: 

^ L  =  Us  cos  C  -  qs  sin  C  (30) 


and 


Drs 

Dt 


=  us  sin  C  +  7s  cos  £ 


(31) 


where  z* ,  rs  are  coordinates  of  nodes  lying  on  the  interface.  These  equations  are  integrated  with  a  4-th  order 
Runge-Kutta  integration  scheme  to  solve  for  the  surface  evolution  in  time.  The  boundary  layer  is  tracked 
with  respect  to  nodes  on  the  jet  surface,  so  that  the  location  of  the  inviscid  core  region  is  given  by: 


*c  =  zs+ 5  sin  C 

rc  =  rs  —  S  cos  C  (32) 


where  (zc,rc)  are  the  coordinates  of  the  node  at  the  edge  of  the  boundary  layer,  and  C  is  the  local  slope  of 
the  free  surface  of  the  jet.  The  unsteady  Bernoulli  equation  requires  a  modification  since  the  nodes  at  the 
edge  of  the  inviscid  core  region  are  no  longer  moving  with  their  instantaneous  velocity.  Neglecting  gravity, 
Bernoulli’s  equation16  can  be  written  as: 


dt  2V 


K 

We 


(33) 


where  <j>  is  the  velocity  potential,  k  is  the  surface  curvature,  and  We  is  the  Weber  number  as  defined  in  the 
Nomenclature.  Performing  an  Eulerian-Lagrangian  transformation  for  nodes  moving  with  the  local  surface 
velocity  gives: 

£->«■+-(£)+*(£)-&  <34> 


where 


and 


d<f)  d<f>  ~  A 

—  =  —  cos/?-gsin/? 
dz  ds 

dd>  dd>  .  „  _ 

_  =  _sm/?+(?coS/? 


(35) 


(36) 


with  /?  the  local  slope  of  the  boundary  layer  interface  with  respect  to  the  centerline  of  the  jet.  In  Equation  34, 
uz  and  ur  are  the  nodal  velocities  on  the  free  surface.  They  are  calculated  by  the  first-order  forward 
differences: 

n  -u  l  n. 

(37) 


,.n  +  l  _ 


rn  +  l 


At 


and 


rn+l 
n-f  1  ;  c 


tC  = 


At 


(38) 


where  superscripts  denote  the  time  level.  Surface  curvature  in  Eq.  34  is  calculated  using  fourth-order 
centered  differences16. 
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Solution  Procedure 

We  use  the  following  procedure  to  solve  for  the  evolution  of  a  finite-length  liquid  jet  and  the  subsequent 
development  of  the  boundary  layer  inside  the  orifice  and  along  the  surface  of  the  jet: 

•  Update  the  integrated  quantities  us,zs ,  and  rs  on  the  jet  surface  and  <j>  at  the  inner  edge  of  the 
boundary  layer  using  a  Runge-Kutta  integration  scheme  on  Equations  25,  30,  31,  and  34; 

•  At  every  node  along  the  free  surface,  update  quantities  dependent  on  the  surface  geometry,  includ¬ 
ing  the  slopes  (3  and  £,  the  curvature  ac,  the  pressure  along  the  free  surface  p,  and  the  derivatives 
dp/dx,  dro/dx  and  dC/dx] 

•  Solve  for  the  inviscid  flow  in  the  core  region  using  the  BEM  solver  for  Laplace’s  equation18.  This 
procedure  returns  values  of  the  velocity  potential  <j> ,  and  normal  velocity  g,  at  the  boundary  of  the 
entire  inviscid  region; 

•  Calculate  properties  at  the  edge  of  the  boundary  layer,  including  the  velocities  d<j>/ds ,  ue  and  d<j)/dy, 
the  pressure  gradient  parameter  A,  and  the  derivatives  dq/ds ,  due/ds ,  dus/ds ,  dv/dx,  dqs/dx,  dp/dx 
and  dk/dx ; 

•  Solve  for  the  boundary  layer  thickness  inside  the  orifice  using  Equation  12  by  marching  from  the  orifice 
inlet  using  a  specified  value  of  S  at  the  inflow  boundary; 

•  Solve  for  the  boundary  layer  thickness  along  the  free  surface  of  the  jet  using  Equation  21  by  marching 
from  the  orifice  exit  to  the  jet  tip; 

•  Calculate  dS fdx  with  a  5-point  difference  scheme  on  the  calculated  values  of  S ; 

•  Solve  for  the  surface  normal  velocity  using  Equation  27;  and 

•  Update  the  location  of  the  edge  of  the  inviscid  core  region  using  Equation  32. 

Initial  Conditions 

Since  the  transition  from  a  no-slip  to  a  free  surface  boundary  condition  is  an  important  element  of  this  model, 
the  inflow  boundary  is  placed  at  a  specified  distance  into  the  nozzle.  The  computational  domain  is  shown  in 
Figure  2,  with  nodes  at  the  interior  edge  of  the  boundary  layer  omitted  for  clarity.  While  this  figure  depicts 
a  simple  cylindrical  orifice  interior,  arbitrary  axisymmetric  contours  can  be  considered  with  the  model.  The 
boundary  layer  thickness  at  the  entry  to  the  domain  ( S0 )  is  a  required  input  for  the  analysis. 

Since  the  choice  of  initial  conditions  essentially  determines  the  nature  of  the  results,  starting  the  calcu¬ 
lation  is  a  key  part  of  this  model.  A  natural  choice  is  to  assume  that  the  jet  is  at  rest  at  t  =  0,  and  turn 
on  the  inflow  at  t  =  0+  with  some  specified  function.  A  problem  with  this  method,  however,  is  that  if  any 
fluid  exists  outside  the  orifice  before  the  calculation  begins,  the  fluid  exiting  the  orifice  tends  to  create  a 
bulging  on  the  jet  surface  just  outside  the  orifice  (klystron  effect).  Although  the  formation  of  surface  waves 
just  after  the  orifice  is  a  desired  result  for  this  model,  we  are  not  interested  in  waves  formed  by  the  klystron 
effect,  but  rather  in  waves  due  to  the  presence  of  the  boundary  layer.  If  the  jet  begins  with  very  little  fluid 
outside  the  orifice,  there  may  be  a  large  computational  cost  involved  in  marching  the  solution  to  a  “steady 
state,”  as  well  as  problems  in  calculating  the  boundary  layer  thickness  on  the  free  surface  of  the  jet. 

The  calculations  performed  in  previous  inviscid  models16  began  with  a  cylindrical  length  of  fluid  with  a 
hemisphere  endcap  moving  at  a  uniform  axial  velocity.  If  the  viscous  calculations  start  with  similar  initial 
conditions,  the  initial  solution  of  the  boundary  layer  thickness  in  the  nozzle  gives  a  step  change  in  profile 
which  leads  to  numerical  instabilities  in  solving  for  the  boundary  layer  thickness  on  the  free  surface.  One 
way  to  alleviate  this  problem  is  to  begin  the  calculation  with  Re  «  oo  and  S  =  0  at  the  inflow  boundary. 
After  beginning  the  calculation,  the  Reynolds  number  is  allowed  to  decrease  to  the  desired  value  over  some 
period  of  time,  and  the  boundary  layer  thickness  at  the  inflow,  tf0,  allowed  to  increase  to  the  desired  value 
over  the  same  period. 

Physically,  this  starting  procedure  amounts  to  a  “ramping”  of  the  viscosity  from  zero  (inviscid  case)  up  to 
the  desired  value  over  a  finite  startup  interval.  While  surface  waves  are  generated  as  a  result  of  this  starting 
transient,  the  waves  are  convected  downstream  such  that  the  near-orifice  behavior  at  t  »  0  is  independent 
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of  the  particular  starting  transient  selected  for  the  simulation.  During  the  startup  interval,  defined  as  A tp, 
So  is  assumed  to  vary  linearly,  and  the  Reynolds  number  is  assumed  to  decrease  according  to  the  following 
function: 

log10  Re  =  ( }aTJ  logl°  Rei  +  \  AiJ  1oSi°  Re°  (39) 

where  A =  ti  -  t0l  and  Re0  and  Rei  are  the  initial  and  final  values  of  the  Reynolds  number,  respectively. 
At  time  tf0,  the  surface  velocity  in  the  nozzle  transitions  from  us  =  1  to  us  =  0  instantaneously,  and  the 
initial  Reynolds  number  is  chosen  as  Reo  =  1012. 


Results 

A  series  of  calculations  were  performed  to  insure  that  the  near-orifice  exit  behavior  was  insensitive  to  the 
assumed  initial  start  transient  at  times  much  greater  than  the  start-up  interval.  Therefore,  the  viscosity¬ 
ramping  technique  was  employed  for  all  results  presented  here.  In  addition,  a  series  of  simulations  were 
conducted  to  verify  that  results  were  insensitive  to  both  temporal  and  spatial  discretization  step  sizes.  A 
typical  calculation  presumes  a  nozzle  internal  length  of  1  jet  radius,  an  initial  jet  length  of  4  radii,  and 
uses  a  time  step  of  At  =  0.002.  Algebraic  stretching16  is  used  to  place  the  15  nodes  in  the  nozzle  with 
points  concentrated  more  heavily  near  the  orifice  exit  plane.  The  grid  spacing  on  the  jet  is  As  =  0.05  after 
the  orifice  for  2  jet  radii,  and  then  transitions  to  As  =  0.10.  Grid  function  convergence  studies  have  been 
conducted  to  insure  that  results  are  insensitive  to  both  temporal  and  spatial  step  sizes  thereby  demonstrating 
that  solutions  are  independent  of  these  parameters.  Algebraic  stretching  and  two  different  grid  spacings  are 
used  as  a  method  of  reducing  the  total  number  of  nodes  and  therefore  the  computational  time.  Typical  run 
times  were  on  the  order  of  a  few  days  on  an  IBM  RISC/System  6000  machine. 

Figure  3  depicts  the  typical  free-surface  evolution  in  the  region  near  the  orifice  exit  during  and  immediately 
after  the  start-up  interval.  Note  that  the  vertical  scale  is  amplified  dramatically;  wave  slopes  are  quite  small 
when  plotted  on  a  true  scale.  This  calculation  was  performed  with  We  =  104,  Re  =  104  and  So  =  0.05.  The 
viscosity  is  turned  on  at  t  =  0.01  time  units,  with  a  startup  interval  of  A =  2  time  units.  As  the  viscosity 
is  increased  during  the  start  transient,  a  surface  wave  is  formed.  This  disturbance  is  convected  downstream 
at  the  mean  velocity  of  the  jet.  While  the  wave  associated  with  the  start  transient  grows  with  time  and 
convects  downstream,  the  surface  near  the  orifice  exit  takes  on  a  steady  shape  with  a  “swelling”  in  diameter 
as  compared  to  that  of  the  orifice.  Physically,  the  swelling  results  from  the  fact  that  the  growth  of  the 
boundary  layer  inside  the  nozzle  serves  as  a  “vena  contracta” .  Immediately  outside  the  nozzle,  an  expansion 
process  takes  place  in  which  the  edge  of  the  jet  adjusts  to  the  local  pressure  in  the  inviscid  core.  This 
phenomena  is  observed  during  die-casting  manufacturing  operations,  and  is  known  as  die  swell.  Numerous 
researchers27-29  have  investigated  this  problem  numerically  for  steady,  2-D  flows. 

When  the  value  of  A t^  is  changed,  surface  waves  associated  with  the  start  transient  are  affected,  but 
not  the  swelling  exhibited  after  the  start  transient  disturbance  convects  downstream.  The  surface  shape 
and  boundary  layer  profile  is  shown  in  Figure  4  for  A =  1,2  and  4  time  units.  When  A is  decreased, 
the  amplitude  of  the  surface  wave  increases.  Note  that  the  assumed  ramping  generates  a  waveform  with 
two  peaks.  The  peaks  exhibit  a  spacing  of  about  one  orifice  radius;  the  last  peak  formed  has  an  amplitude 
roughly  double  that  of  the  initial  crest.  This  trend  is  apparent  at  very  short  start  intervals  and  ultimately 
leads  to  waves  which  are  too  steep  to  resolve  with  the  assumed  boundary  layer  profiles.  Since  the  waves  form 
near  the  end  of  the  ramp  interval,  the  shorter  intervals  form  the  wave  earlier  in  the  calculation  and  the  wave 
has  longer  to  move  downstream  in  the  figure.  The  boundary  layer  thickness  and  jet  surface  shapes  appear 
nearly  identical  after  the  wave  has  been  convected  downstream  from  the  orifice  as  discussed  previously. 

Figure  5  shows  how  the  surface  and  boundary  layer  edge  velocities  change  in  the  nearfield  region  over 
the  first  5.5  time  units  of  the  calculations.  For  A tp  —  1  time  unit,  the  edge  velocity  exhibits  a  few  ripples 
(associated  with  the  more  violent  start  transient),  but  essentially  remains  constant  at  ue  =  1  along  the  surface 
at  all  times.  The  distance  over  which  the  surface  velocity  transitions  from  the  no-slip  condition  to  us  =  1 
grows  steadily  with  time,  and  is  shown  to  be  the  same  in  all  three  cases.  Very  near  the  orifice  exit  (say,  for 
z  <  1),  results  begin  to  asymptotically  approach  a  fixed  us  value.  Unfortunately,  computational  restrictions 
prohibit  simulations  to  greater  times  in  which  the  entire  nearfield  flow  becomes  steady.  Fortunately,  the 
free  surface  development  occurs  in  the  region  very  near  the  orifice  exit  (z  <  0.5),  such  that  the  current 
calculations  reveal  at  steady  behavior  in  this  region  as  highlighted  in  Fig.  4  and  subsequent  results. 
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The  growth  of  the  boundary  layer  thickness  for  a  A =  2  time  unit  ramp  interval  is  shown  in  Figure  6. 
The  boundary  layer  grows  rapidly  during  the  start  transient,  but  after  about  3  time  units,  the  behavior  is 
invariant  with  time.  In  addition,  note  the  dramatic  thinning  of  the  boundary  layer  in  this  region  immediately 
downstream  of  the  orifice  exit.  At  a  distance  of  only  1/2  a  jet  radius  downstream  of  the  orifice,  the  boundary 
layer  has  all  but  vanished.  Keeping  in  mind  that  a  substantial  velocity  difference  exists  between  surface 
and  edge  velocities  at  this  spatial  location,  the  boundary  layer  is  actually  transitioning  to  a  shear  layer  in 
the  region  very  near  the  orifice  exit.  Viscous  forces  would  presumably  lead  to  a  broadening  of  the  shear 
layer  as  the  flow  progresses  downstream.  Since  the  Reynolds  numbers  in  the  simulations  are  very  large,  the 
relaxation  lengths  are  also  very  large  and  cannot  be  resolved  due  to  computational  time  restrictions. 

Figure  7  shows  the  actual  boundary  layer  profiles  near  the  orifice  exit  once  steady-state  conditions  have 
been  attained  ( t  «  5.5).  Once  again,  note  that  the  vertical  scale  has  been  amplified  dramatically  for 
presentation  purposes.  In  the  upper  plot,  the  abcissa  represents  the  velocity  relative  to  the  local  surface 
velocity  and  the  ordinate  is  the  distance  from  the  surface  (essentially  the  —r  direction).  As  expected,  the 
boundary  layer  grows  until  the  orifice  exit  (Point  “3”  in  the  figure).  Outside  the  orifice,  the  rapid  thinning  of 
the  boundary  layer  is  evident,  and  by  Point  “7”  the  layer  is  very  thin.  The  thinning  of  the  boundary  layer  to 
a  shear  layer  (vortex  tube),  is  known  to  be  an  inherently  unstable  phenomenon.  As  the  shear  layer  develops, 
the  presumed  boundary  layer  velocity  profile  (Eq.  14)  becomes  erroneous  in  that  it  cannot  properly  model 
rotation  within  the  boundary  layer.  Therefore,  there  will  be  a  tendency  toward  instability  in  the  region 
far  downstream  of  the  orifice  exit  plane.  In  spite  of  this  fact,  the  near  field  solution  is  resolved  adequately 
using  Eq.  14  velocity  profiles.  More  research  is  required  to  determine  whether  or  not  viscous  dissipation  will 
ultimately  limit  the  formation  of  this  unstable  vortex  tube. 

Effect  of  So,  Re  and  We 

A  series  of  calculations  were  performed  to  quantify  the  amount  of  swelling  (as  a  percentage  of  orifice  radius) 
for  various  initial  conditions.  These  calculations  were  performed  at  We  =  104  for  various  Jo  and  Re  values. 
Analysis  of  the  results  indicated  that  the  swelling  correlated  most  closely  with  the  thickness  of  the  boundary 
layer  at  the  exit  of  the  orifice  passage.  Figure  8  presents  a  summary  of  swelling  results  plotted  against 
this  thickness.  Results  indicate  that  the  amount  of  swelling  is  essentially  independent  of  Reynolds  number 
(except  for  its  obvious  implications  in  setting  the  boundary  layer  thickness  at  the  exit  plane).  Assuming  that 
all  the  swelling  takes  place  within  about  the  first  1/2  jet  radius  (which  is  consistent  with  our  calculations), 
the  net  “cone  angles”  corresponding  to  the  swelling  results  in  Fig.  8  lie  in  the  range  of  1-5°.  These  small 
angles  would  be  difficult  to  observe  experimentally.  In  addition,  the  bulk  of  experiments  with  images  near 
the  orifice  exit  are  for  conditions  where  substantial  gas-phase  pressure  interactions  are  present. 

Simulations  performed  at  Reynolds  numbers  other  than  104  show  similar  behavior  to  that  noted  in  Figs. 
3-7.  In  general,  the  waves  associated  with  the  start  transient  become  more  prominent  at  lower  Re  values, 
but  the  steady-state  behavior  is  quite  insensitive  to  this  parameter  (for  a  given  value  of  S  at  the  orifice  exit) . 
Results  indicated  an  insensitivity  to  Weber  number  as  well;  nearly  identical  surface  shapes  were  obtained  for 
a  range  of  Weber  numbers  for  simulations  employing  fixed  J0  and  Re  values.  For  these  reasons,  the  swelling 
phenomena  can  effectively  be  characterized  by  So  alone;  the  amount  of  viscosity  present  is  only  important 
in  as  much  as  it  affects  this  parameter. 


Discussion 

Assessing  the  results  shown  in  the  previous  section  in  the  light  of  experimental  observations  of  other  re¬ 
searchers  as  described  in  the  introduction,  we  can  draw  several  important  conclusions.  First  we  note  that 
the  present  study  has  shown  that  steady-state  solutions  do  exist  for  a  high-speed  laminar  jet  in  the  absence 
of  pressure  perturbations  from  the  gas  phase.  This  result  indicates  that  atomization  mechanisms  which  rely 
on  boundary  layer  relaxation  or  velocity  profile  rearrangement  cannot  generate  the  instabilities  required  to 
atomize  the  jet.  Here,  we  note  that  the  experiments  of  Reitz  and  Bracco3  and  others  also  support  this 
conclusion. 

The  swelling  phenomena  depicted  in  all  calculations  with  this  model  can  aid  in  explaining  many  of  the 
experimental  observations.  In  the  presence  of  a  gas-phase,  the  swelling  of  the  jet  provides  a  large-amplitude 
disturbance  for  initiation  of  aerodynamic  instabilities  which  are  a  critical  element  of  Taylor’s  model.  For  the 
high-speed  jet,  these  instabilities  could  be  amplified  so  fast  that  atomization  occurs  in  the  region  very  near 
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the  exit  plane.  Disturbances  from  these  atomization  events  could  be  fed  back  to  the  exit  plane  such  that  a 
steady-state  solution  no  longer  exists.  This  scenario  is  quite  plausible  for  a  realistic  situation  in  which  the 
“start  transient”  is  more  chaotic  and  three-dimensional  in  nature  than  the  idealized  approach  used  in  our 
modeling. 

Since  the  amount  of  swelling  correlates  well  with  the  boundary  layer  thickness  at  the  exit  plane,  we  can 
also  draw  some  conclusions  about  the  effect  of  boundary  layer  turbulence  and  orifice  design.  Inasmuch  as 
turbulence  reduces  the  thickness  of  the  boundary  layer  (as  compared  to  laminar  case),  the  amount  of  jet 
swelling  would  tend  to  be  reduced,  and  the  corresponding  perturbations  from  the  gas-phase  would  also  be 
smaller.  On  this  basis,  the  high-speed  jet  (where  gas-phase  pressure  distribution  has  a  major  effect)  would 
be  more  stable  under  the  case  of  a  turbulent  boundary  layer.  This  conclusion  agrees  with  experimental 
observations  of  atomizing  jets3,9. 

In  the  low-speed  regime,  turbulence  could  play  a  significant  role.  Since  our  solutions  (Figs.  3-8)  presume 
no  gas-phase  interaction  and  a  laminar  boundary  layer,  we  predict  that  steady  solutions  could  be  possible 
in  the  low-speed  regime  (particularly  for  orifice  designs  generating  thin  boundary  layers  at  the  exit).  Under 
these  conditions,  forces  due  to  turbulent  fluctuations  in  the  radial  direction  could  be  the  dominant  mechanism 
destabilizing  the  jet.  These  observations  are  consistent  with  the  experimental  measurements7,8  which  were 
conducted  in  a  lower-speed  regime. 

As  the  orifice  passage  is  lengthened  (with  fixed  flowrate) ,  the  boundary  layer  at  the  exit  plane  thickens 
and  the  swelling  phenomena  becomes  more  pronounced.  On  this  basis,  one  would  predict  that  increasing  Lj D 
would  cause  the  jet  to  be  more  unstable  -  a  result  consistent  with  experimental  observations  of  McCarthy  and 
Molloy4  for  the  case  of  a  low-speed  jet.  However,  for  high-speed  jets,  Reitz  and  Bracco  observed  the  opposite 
trend  with  increased  L/D.  As  pointed  out  by  these  (and  other)  authors,  cavitation  may  play  a  substantial 
role  complicating  the  behavior  since  it  can  provide  a  mechanism  for  generating  freestream  turbulence6  which 
is  not  included  in  our  model.  In  addition,  very  long  nozzles  tend  to  damp  out  flow  perturbations  caused 
by  surface  imperfections  at  the  orifice  inlet.  Recent  research31  has  shown  this  phenomena  to  be  important 
in  determining  jet  atomization  characteristics  in  some  cases.  Obviously,  this  mechanism  could  also  obscure 
conclusions  made  regarding  orifice  L/D  effects. 


Conclusions 

A  model  has  been  developed  to  describe  the  nonlinear,  unsteady  evolution  of  an  axisymmetric  liquid  jet 
which  contains  a  thin  boundary  layer.  Results  from  simulations  indicate  that  steady-state  solutions  do  exist 
for  the  case  of  a  laminar  boundary  layer  with  negligible  gas-phase  interactions.  Under  these  conditions, 
the  jet  swells  a  few  percent  in  diameter  and  the  boundary  layer  thins  to  a  shear  layer  within  the  first  1/2 
orifice  radius  of  the  exit  plane.  These  results  demonstrate  that  atomization  mechanisms  which  rely  on 
boundary  layer  relaxation  cannot  explain  this  process.  Moreover,  the  swelling  phenomena  could  provide  a 
large-amplitude  perturbation  for  amplification  of  gas-phase  pressure  interactions  and  explain  experimental 
observations  in  both  high-speed  (atomizing)  and  low-speed  regimes. 

The  behavior  predicted  by  the  model  is  consistent  with  experimental  observations  regarding  orifice  Lj D 
effects  in  the  case  of  low-speed  flows.  However,  discrepancies  do  exist  when  comparing  with  the  experiments 
on  L/D  effects  in  higher  speed  flows  consistent  with  the  atomization  regime.  These  discrepancies  might  be 
explained  by  cavitation  or  surface  imperfections  which  complicate  this  issue. 
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Nomenclature 

e  =  unit  vector 

p  —  pressure 

q  =  normal  velocity 

r  =  radial  coordinate 

Re  =  Reynolds  number 

s  =  distance  along  the  surface 

t  —  time 

u  —  velocity  parallel  to  the  wall  or  free  surface 

v  =  velocity  perpendicular  to  the  wall  or  free  surface 

We  =  Weber  number,  We  =  pU2a/a 

x  —  coordinate  parallel  to  the  wall  or  free  surface 

y  =  coordinate  perpendicular  to  the  wall  or  free  surface 

z  =  axial  coordinate 

/?  =  slope  of  the  boundary  layer  interface 

S  =  boundary  layer  thickness 

S*  =  displacement  thickness 

rj  —  boundary  layer  parameter,  0  <  rj  =  y/S  <  1 

k  =  surface  curvature 

A  =  pressure  gradient  parameter 

p  —  viscosity 

<j)  —  velocity  potential 

( Tij  =  stress  tensor 

0  —  momentum  thickness 

Subscripts: 

e  =  edge  of  the  boundary  layer 
s  =  surface 
c  =  core  region 
Superscripts: 
n  =  time  level 
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Appendix 

Detailed  equations  for  the  coefficients  of  the  momentum  equations  for  a  solid  wall  (Eq.  12)  and  along  the 
free  surface  (Eq.  21)  are  provided  herein. 


Solid  Wall  Boundary  Layer  Equations 

By  expanding  derivatives  and  integrals  in  Eq.  5,  the  integral  form  of  the  momentum  equation  can  be  written: 

A80  +  u\8 jjL  +  ue8^8*  +  B8 2  —  C83  =  D 


where: 


Next,  we  substitute  for  displacement  and  momentum  thicknesses  (Eqs.  9,10)  and  use  the  backward  difference 
(Eq.  11)  as  a  local  approximation  to  dS/dx.  Following  this  procedure,  we  group  terms  to  obtain  the  form 
shown  in  Eq.  12.  The  resulting  coefficients  can  be  written: 


c0  =  D 


where: 
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Free  Surface  Boundary  Layer  Equations 

In  this  case,  the  original  form  of  the  momentum  equation  can  be  written  (after  substitution  for  momentum 
and  displacement  thicknesses)  as: 


r  f2  as  Sds  ds  , 

a0  +  aiS  +  a2S  +a3—  +  a,S-  +  a5Tt+a6S- 


=  0 


Upon  application  of  a  backward  difference  in  both  space  and  time,  the  equation  can  be  rearranged  to  the 
result  shown  in  Eq.  21.  Here,  the  f>,  values  can  be  expressed: 
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where  the  coefficients  a0  —  a6  are: 
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The  velocity  derivative  dv/dx  can  be  written  in  terms  of  derivatives  of  <j>  and  q ,  as  well  as  local  surface 
geometry  (£,  /3).  Other  spatial  derivatives  appearing  in  these  equations  are  evaluated  using  centered  dif¬ 
ferencing  since  at  a  given  time,  arguments  in  these  derivatives  are  known  along  the  entire  surface.  Time 
derivatives  appearing  in  the  above  relations  are  evaluated  using  a  backward  difference  as  utilized  for  the 
unknown  8 . 
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Figure  Captions 


Figure  1.  Axisymmetric  Boundary  Layer  Nomenclature. 

Figure  2.  Computational  Domain  for  the  Finite-Length  Viscous  Jet.  (Nodes  Along  Inner  Edge  of  Bound¬ 
ary  Layer  are  not  Shown  for  Clarity) 

Figure  3.  Boundary  Layer  Development,  We  =  104,  Re  =  104,  A tp  =  2.0.  Here  the  Upper  Curves  Repre¬ 
sent  the  Jet  Surface  While  the  Lower  Curves  Denote  the  Inner  Edge  of  the  Boundary  Layer. 

Figure  4.  Jet  Profile  Near  the  Orifice  Exit  for  We  =  104,  Re  =  104  and  t  =  5.8  Time  Units. 

Figure  5.  Edge  and  Surface  Velocities,  We  =  104,  Re  =  104,  and  At  =  0.5  Between  Curves,  (a)  A tp  =  1, 
(b)  A tp  =  2,  (c)  At ^  =  4  Time  Units. 

Figure  6.  Boundary  Layer  Thickness,  We  =  104,  Re  =  104,  A tp  =  2.0,  with  A^  =  0.5  Between  Curves. 

Figure  7.  Boundary  Layer  Velocity  Profiles  at  Various  Locations  Along  the  Jet,  We  =  104,  Re  =  104. 

Figure  8.  Jet  Swelling  as  a  Function  of  Boundary  Layer  Thickness  at  the  Orifice  Exit  Plane  and  as  a 
Percentage  of  Orifice  Radius,  We  =  104. 


28 


Radial  Coordinate 


Axial  Coordinate 


Figure  4:  Jet  Profile  Near  the  Orifice  Exit  for  We  =  104,  Re  =  104  and  t  =  5.8  Time  Units. 
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Figure  8:  Jet  Swelling  as  a  Function  of  Boundary  Layer  Thickness  at  Orifice  Exit  Plane  as  a  Percentage  of 
Orifice  Radius,  We  =  104. 


6  Appendix  B  -  F-l  Engine  Simulations 
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CONTRIBUTIONS  OF  ATOMIZATION  TO  F-l 
ENGINE  COMBUSTION  INSTABILITIES 

S.  D.  Heister  *  and  C.  Chao  t 
January  28,  1999 


Nomenclature 

a  =  orifice  radius 
P  =  pressure 

q  =  velocity  normal  to  local  surface 
t  —  time 
U  =  velocity 

We  =  Weber  number,  IFe  =  pgUga/cr 

x  =  horizontal  coordinate 

y  =  vertical  coordinate 

/c  =  surface  curvature 

p  =  density 

e  =  density  ratio,  pg/p 

a  =  surface  tension 

ujg  =  acoustic  oscillation  frequency 

=  natural  frequency  of  column  (Eq.  1) 
(j)  =  velocity  potential 

Subscripts 
()ff  =  gas  phase 


Abstract 

Numerical  simulations  have  been  performed  to  assess  the  nonlinear,  unsteady  behavior  of  fuel  jets  emanating 
from  an  injector  in  the  presence  of  a  transverse  (tangential)  acoustic  wave.  Conditions  for  the  simulations 
have  been  set  to  pertain  to  injector  designs  used  during  the  development  of  the  F-l  liquid  rocket  engine. 
Results  indicate  that  injector  designs  which  were  unstable  have  substantial  coupling  between  the  acoustic 
wave  and  the  natural  frequency  of  oscillation  of  the  fuel  column  emanating  from  the  orifice.  This  instability 
scenario  is  discussed  in  the  context  of  other  theories  offered  as  explanation  for  F-l  instabilities. 
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Introduction 


During  the  development  of  the  F-l  liquid  rocket  engine  which  was  used  as  first  stage  propulsion  on  the 
Saturn  V  launch  vehicle,  tangential  mode  instabilities  were  observed  in  a  variety  of  injector  designs.  Literally 
thousands  of  tests  were  conducted1,  providing  the  most  extensive  fullscale  database  available  even  today. 
For  this  reason,  the  F-l  engine  continues  to  be  a  subject  of  study;  providing  a  substantial  amount  of  engine 
response  data  with  which  to  correlate  new  experiments  and  analytic/numerical  models. 

The  F-l  engine  utilized  liquid  oxygen  and  RP-1  (kerosene)  propellants  in  an  impinging- type  injector 
design  depicted  schematically  in  Fig.  1.  Fuel  was  injected  in  an  impinging  doublet  pattern  in  which  the 
angle  between  the  jets  was  typically  40° ,  while  the  LOX  was  injected  in  a  triplet  pattern  in  which  an  axial 
center  jet  is  mixed  with  two  canted  jets  at  the  impingement  point.  The  angle  between  the  canted  jets  was 
typically  kept  at  40°  as  in  the  case  of  the  fuel  jets.  Various  injector  designs  were  tested  in  the  program;  main 
differences  between  injectors  include  the  orifice  size/number  and  the  presence  (or  absence)  of  baffles.  Many 
of  the  engines  using  the  early  injector  designs  exhibited  instabilities  in  the  450-550  Hz  range  consistent  with 
the  first  tangential  mode  for  this  large  engine. 

In  recent  years,  several  theories  have  been  adopted  to  explain  the  presence  of  these  instabilities  in  F- 
1  and  other  impinging  element  injector  designs.  Most  of  these  theories  focus  on  the  behavior  of  the  fuel 
jets  because  they  are  the  least  volatile  of  the  two  propellants  and  most  likely  to  provide  a  rate-controlling 
mechanism  for  combustion  processes.  A  stability  correlation  after  Hewitt2  indicates  that  the  maximum 
sustainable  instability  frequency  decreases  as  the  injection  parameter  a/Ui  is  increased.  Here,  U{  is  the  fuel 
injection  velocity  and  a  is  the  corresponding  orifice  radius.  The  basic  jet  atomization  frequency  resulting 
from  the  impingement  process  has  been  shown  to  be  of  the  same  order  as  those  observed  in  testing2,3,  thereby 
providing  one  mechanism  in  agreement  with  Hewitt’s  criteria.  In  addition,  flame  straining  and  extinction 
times  have  been  shown4  to  provide  frequencies  in  accordance  with  the  observed  trends  and  those  predicted 
by  Hewitt’s  criteria. 

Finally,  fuel  jet  aerodynamic  excitation5  has  been  suggested  as  a  potential  contributor  since  the  column 
of  liquid  has  a  natural  frequency  of  excitation  which  may  be  near  that  of  the  acoustic  frequency.  In  this  paper, 
we  investigate  this  latter  theory  with  some  detailed  simulations  for  both  stable  and  unstable  F-l  injector 
designs.  The  next  section  provides  a  brief  description  of  the  model,  followed  by  results  and  discussion  in  the 
context  of  the  other  theories  described  above. 


Model  Description 

The  two-dimensional  model  used  in  the  simulations  conducted  herein  makes  use  of  Boundary  Element  Meth¬ 
ods  (BEMs)  to  provide  high-resolution,  time-accurate  modeling  of  a  liquid  column  subjected  to  acoustic 
forces  from  the  gaseous  phase.  Since  we  have  described  details  of  the  fully-coupled,  nonlinear  models  in 
previous  work5-8,  we  shall  provide  only  a  brief  discussion  in  the  present  study.  The  model  used  for  the 
present  simulations  presumes  incompressible,  inviscid  flow  in  both  liquid  and  gaseous  phases.  Since  the 
orifice  diameter  is  typically  2-3  orders  of  magnitude  less  than  the  chamber  diameter,  an  acoustic  wave  can 
effectively  be  represented  as  a  time- varying  incompressible  flow;  i.e.  spatial  variations  across  the  tiny  jet 
length  scale  are  assumed  to  be  negligible.  The  inviscid  assumption  is  necessary  for  use  of  BEMs;  transient 
flow  separation  present  in  a  complete  viscous  solution  introduces  formidable  complications  in  the  modeling. 

We  choose  the  jet  radius  (a),  liquid  density  (p),  and  a  characteristic  peak  acoustic  velocity  ( Ug )  as 
dimensions.  Under  this  nondimensionalization,  the  Weber  number,  We  =  pgUga/cr  and  the  gas/liquid 
density  ratio,  e  =  pgj p  are  two  dimensionless  variables  characterizing  these  flows.  A  third  parameter  which 
is  important  in  these  unsteady  simulations  is  the  natural  (lowest  order)  excitation  frequency  of  the  column. 
In  Ref.  5,  a  linear  analysis,  similar  to  that  conducted  by  Lamb  for  vibrations  of  a  droplet,  provides  the 
natural  frequency  of  the  fuel  jet  under  transverse  excitation: 
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The  ratio  of  the  imposed  acoustic  frequency  (u>g)  and  u)n  represents  another  important  parameter  in  this 
particular  problem. 
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For  an  inviscid  gas  or  liquid  domain,  velocity  potentials  0ff ,  0  exist  and  must  satisfy  Laplace’s  equation: 

V20  =  V20ff  =  0  (2) 


The  unsteady  Bernoulli  equation  provides  conditions  relating  velocity  potentials  at  the  gas/liquid  interface: 


^  +  i(V0)2  +  Ps  +  ^-O 


(3) 


where  Pg  is  the  dimensionless  gas  pressure  at  the  interface,  and  k  is  the  local  surface  curvature.  On  the  gas 
side  of  the  interface,  Bernoulli’s  equation  is: 


^  +  |(v^)2  +  JP</  =  °  (4) 

Mathematically,  Eqs.  2-4  provide  a  well-posed  set  of  equations  for  velocity  potentials,  the  gas  pressure 
at  the  interface,  and  the  shape  of  the  interface  (implicit  in  the  surface  curvature,  k).  This  set  of  equations 
is  solved  using  a  BEM  which  begins  with  an  integral  representation  of  Laplace’s  equation: 

a<Kfi)  +  -  qG]dT  =  0  (5) 


where  0(rj)  is  the  value  of  the  potential  at  a  point  rj,  T  denotes  the  boundary  of  the  domain,  and  G  is  the 
free-space  Green’s  function  corresponding  to  Laplace’s  equation.  An  analogous  form  of  Eq.  5  can  also  be 
derived  for  the  gas  phase  potential.  For  a  well-posed  problem,  either  0  or  q  —  d<t>/dn  must  be  specified  at 
each  “node”  on  the  boundary.  Here  n  is  the  outward  normal  to  the  boundary  so  that  q  represents  the  velocity 
normal  to  the  boundary.  The  quantity  a  in  Eq.  5  results  from  singularities  introduced  as  the  integration 
passes  over  the  boundary  point,  r*. 

In  the  case  of  a  2-D  flow  (letting  x  and  y  represent  the  coordinates),  the  free-space  Green’s  function  may 
be  written: 

G=  —  ln|r-rf|  =  ^  ln[(x  -  x,)2  +  (y  -  j/i)2]  (6) 

We  presume  that  both  0  and  q  vary  linearly  along  the  length  of  a  given  element.  This  assumption  permits 
the  construction  of  a  set  of  matrices  involving  the  nodal  values  of  0  and  q  and  the  integrals  (over  a  given 
element)  given  in  Eq.  5.  For  this  linear  element  case,  the  integration  across  a  segment  can  be  carried  out 
analytically.  Singularities  resulting  from  integration  across  a  segment  containing  the  base  point  are  also 
integrable6. 

The  main  challenge  in  developing  models  capable  of  tracking  large  deformations  of  an  interface  lies  in 
the  treatment  of  the  free  surface  itself.  Here,  the  use  of  a  BEM  is  desirable  since  nodes  can  be  tracked 
such  that  they  always  lie  on  the  interface;  no  interpolations  of  surface  location  (and  the  inherent  numerical 
errors)  are  required.  Since  capillary  forces  are  important,  it  is  crucial  to  develop  a  treatment  capable  of 
accurately  determining  surface  curvature  at  all  times  during  the  simulation.  For  this  reason,  all  models 
employ  fourth-order  centered  differencing  (on  a  generalized,  variable  spacing  mesh)  to  determine  surface 
curvature7. 

We  track  surface  nodes  along  the  local  liquid  velocity  vector.  Under  this  assumption,  flow  kinematics 
require: 

Dx  _  50  Dy~  _  50 

Dt  ~  dx  Dt  ~~  dy  U 

where  the  notation  D()/Dt  implies  a  Lagrangian  derivative  for  points  on  the  surface  moving  with  the  local 
liquid  velocity.  Recognizing  that  our  BEM  solver  will  return  velocities  normal  to  the  surface,  we  employ  the 
velocity  transformations: 

=  |f sin(P)  +  <1C0S(P)  =  J-SC0SW  ~  qsin{/3)  (8) 

where  j3  is  the  local  wave  slope  and  50/5s  is  the  velocity  tangential  to  the  local  surface.  This  tangential 
velocity  is  calculated  using  5-point  centered  differences  on  surface  values  of  0.  The  local  wave  slope,  /?,  is 
calculated  by  determining  the  slope  of  a  parabola  fit  through  three  points  in  the  region  of  interest9 . 
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Dynamics  of  the  interface  are  addressed  through  the  unsteady  Bernoulli  equation.  The  relations  in 
Eqs.  3  and  4  are  valid  for  an  Eulerian  system  in  which  the  grid  remains  fixed.  Since  we  desire  nodes  on 
the  interface  to  move  in  a  Lagrangian  sense,  an  Eulerian  -  Lagrangian  transformation  is  required.  Here, 
we  “track”  nodes  on  the  interface  by  assuming  that  they  travel  with  the  local  liquid  velocity.  Under  this 
assumption,  the  Eulerian  -  Lagrangian  transformation  can  be  written: 

5L>  =  *Q+W.V()  (9) 

and  the  applicable  forms  of  the  Bernoulli  equations  (Eqs.  3,4)  become: 


D<t> 

Dt 


(10) 


CDt  =  “  P *  (n) 

Equations  7,  10,  and  11  are  integrated  in  time  using  a  fourth-order  Runge-Kutta  scheme.  A  stable, 
consistent  procedure6  has  been  developed  to  handle  the  coupled,  nonlinear  boundary  conditions  at  the 
interface  (Eqns  10,11).  Equation  10  is  integrated  to  update  <j>  values  on  the  interface  which  permits  solution 
of  Laplace’s  equation  in  the  liquid  phase  to  obtain  surface  velocities.  Since  we  tract  nodes  along  this  velocity 
vector,  we  must  have  qg  =  —q  for  gas  nodes  on  the  free  surface.  Using  this  condition,  Laplace’s  equation 
is  solved  in  the  gas  phase  to  determine  <f>g  on  the  surface.  This  current  <j>g  value  is  used  in  a  first  order 
backward  difference  to  estimate  D<j)g/Dt  such  that  Eq.  11  can  be  used  to  solve  to  the  updated  Pg  values. 

Since  the  nodes  on  the  interface  are  allowed  to  move  with  their  local  velocity,  over  time  they  tend  to 
group  themselves  in  regions  of  high  curvature.  This  phenomena  leaves  regions  of  lower  curvature  poorly 
defined.  To  alleviate  this  problem,  the  surface  mesh  is  regridded  using  a  series  of  cubic  splines  (for  surface 
coordinates  and  <f>)  at  each  time  step  to  keep  the  spacing  between  the  nodes  constant  along  the  surface. 


Computational  Grid  and  Boundary  Conditions 

Figure  2  highlights  the  computational  grid  and  boundary  conditions  for  the  simulation  of  jet  behavior  in  the 
presence  of  a  transverse  acoustic  excitation.  The  2-D  assumption  restricts  us  to  considering  a  “slice”  of  the 
jet  upstream  of  the  impingement  point  shown  in  Fig.  1.  The  Weber  number  (We),  gas/liquid  density  ratio 
(e),  and  gas/liquid  frequency  ratio  (u)g/i on)  are  input  parameters  for  a  given  simulation.  We  presume  that 
the  acoustic  oscillation  can  be  represented  by  a  simple  cosine  wave,  so  that  the  gas-phase  velocity  potential 
far  from  the  column  may  be  written: 

(j)g  =  xcos(wgt/2)  (12) 

The  factor  of  1/2  is  included  inside  the  cosine  function  to  account  for  the  fact  that  the  inviscid  solution  is 
insensitive  to  the  direction  of  gas  flow,  i.e.  the  column  will  broaden  under  the  presence  of  an  acoustic  wave 
moving  either  direction.  Numerical  experiments5  using  the  known  analytic  solution  for  uniform  flow  over  a 
cylinder  indicate  that  farfield  conditions  may  be  accurately  assumed  if  the  outer  gas  boundary  is  placed  15 
jet  radii  from  the  origin. 

The  column  is  subjected  to  a  virtual  mass  force  due  to  the  accelerating  external  flow.  Since  there  is 
a  freestream  pressure  gradient  present,  the  flow  is  not  symmetric  about  the  y  axis.  However,  for  low  gas 
density  flows,  similar  analyses  on  droplets10  indicate  that  the  apparent  mass  force  is  negligible.  In  addition, 
we  have  also  run  numerous  simulations  for  the  complete  upper  half  plane5  (as  shown  in  Fig.  2);  at  the 
conditions  of  interest,  the  quarter-plane  results  are  equivalent.  The  advantage  in  using  the  quarter-plane  is 
obvious;  it  halves  the  number  of  nodes  in  the  simulation.  Since  the  matrix  inversion  time  scales  with  the 
cube  of  the  number  of  nodes,  the  quarter-plane  domain  takes  about  1/8  the  time  to  run  as  compared  to  the 
half-plane  domain.  For  this  reason,  we  utilized  the  domain  shown  in  Fig.  2  for  all  simulations  presented 
herein. 

For  liquid  nodes  along  the  assumed  symmetry  axes  in  Figure  2  the  normal  velocity,  q ,  is  assumed  to  be 
zero.  For  the  gas  domain,  symmetry  dictates  that  qg  =  0  along  the  line  y  =  0,  the  assumed  symmetry  along 
the  x  =  0  line  dictates  that  (j)g  be  a  constant  equivalent  to  the  value  prescribed  by  the  freestream  potential 
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Table  1:  F-l  Injector  Design  Data:  Fuel  Orifice  Test  Data 


Injector  Desigation 

Dimm.  mm 

Ug,  Hz 

A  P/Pc  % 

0Jg/un 

5U- Flat  face 

3.66 

538 

150 

3.0 

5U-Baffled 

4.04 

460 

65 

2.9 

Double  Row  Cluster  (DRC) 

2.79 

454 

400 

1.7 

Prelim.  Flight  Rating  Tests  (PFRT) 

5.79 

- 

0 

5.1* 

Flight  Rating  Test  (FRT) 

7.14 

- 

0 

7.1* 

*  Assumes  u>9  =  460  Hz  IT  chamber  mode. 


in  Eq.  12.  The  Bernoulli  equations  (Eqs.  11,12)  prescribe  the  necessary  conditions  for  the  interface.  A 
dimensionless  time  step  of  0.005  was  used  in  all  simulations. 

Typical  grids  employ  21  nodes  on  outer  boundary  of  gas  and  15  gas  nodes  along  each  symmetry  plane 
for  the  gas  phase.  In  the  liquid  phase,  26  nodes  had  been  used  on  the  interface  and  nine  nodes  along  each 
symmetry  plane.  Figure  3  depicts  the  column  aspect  ratio  as  a  function  of  time  for  typical  We  and  e  values 
for  three  different  computational  grids.  In  this  figure,  the  fine  and  coarse  grids  employ  double  and  half  the 
number  of  nodes  respectively,  as  compared  to  the  “basis  grid”  described  above.  The  simulation  demonstrates 
that  the  solution  is  independent  of  grid  provided  we  use  at  least  as  many  nodes  as  that  employed  in  the 
basis  grid.  Using  this  grid,  the  model  can  also  replicate  analytic  solutions  for  uniform  flow  over  a  cylinder 
to  within  1%.  For  these  reasons,  this  grid  has  been  utilized  for  all  results  in  this  work. 

Table  1  highlights  data1  from  both  unstable  and  stable  injectors  tested  during  F-l  engine  development. 
The  Flight  Rating  Test  (FRT)  design  is  essentially  the  configuration  which  was  used  on  actual  flight  units. 
The  most  unstable  injector  tested  was  the  Double  Row  Cluster  (DRC)  which  exhibited  chamber  pressure 
oscillations  as  high  as  400%  of  the  mean  value.  Using  a  dimensional  form  of  Eq.  1,  column  natural  frequencies 
were  calculated  for  each  fuel  orifice  design  uln  =  yJbpgV /[{p  +  Pg)° 3]  using  a  fuel  specific  gravity  of  0.81  and 
surface  tension  of  0.0019  lbf/f.  Using  this  information,  the  frequency  ratio  u )g/un  was  then  computed.  Gas 
density  was  estimated  for  nominal  combustion  conditions  (chamber  pressure  of  1125  psi,  mixture  ratio  of 
2.4)  using  one-dimensional  chemical  equilibrium  calculations.  The  frequency  ratio  for  the  stable  injectors 
was  estimated  assuming  a  first  tangential  mode  frequency  of  460  Hz  for  this  large  chamber. 


Modeling  Results 

To  highlight  differences  between  injectors,  we  decided  to  model  the  most  unstable  (DRC)  and  stable  (FRT) 
injector  designs  in  our  simulations.  Using  nominal  fuel  and  gas  densities,  the  applicable  density  ratio  for 
the  F-l  simulations  is  e  =  0.0064.  Prescribing  an  acoustic  intensity  (in  terms  of  a  Weber  number)  is 
difficult  because  the  exact  acoustic  environment  to  which  the  fuel  jets  are  exposed  is  largely  unknown 
and  not  measured  in  typical  testing.  For  this  reason,  a  baseline  Weber  number  of  0.1  was  assumed  for 
most  simulations.  This  value  corresponds  to  peak  acoustic  velocities  in  the  range  from  1-2  f/s;  modest 
perturbations  which  could  be  consistent  with  early  stages  of  instability. 

In  addition  to  actual  column  shapes,  we  have  found  that  tracking  the  motion  of  the  top  node  ( x  =  0)  on 
the  interface  to  provide  insight  into  the  nonlinear  deformation  processes  depicted  by  the  model.  In  addition, 
the  column  aspect  ratio  (AR),  defined  as  the  ratio  of  column  height  to  its  width,  is  one  variable  which 
quantifies  the  level  of  deformation  of  the  jet.  Prolate  and  oblate  shapes  are  formed  when  AR  >  1  and 
AR  <  1,  respectively. 

Using  the  frequency  ratios  provided  in  Table  1 ,  Figures  4  and  5  depict  the  response  of  DRC  and  FRT 
injector  fuel  jets  in  terms  of  top  node  radial  motion  and  column  aspect  ratio.  In  Fig.  4,  we  see  that  both 
designs  exhibit  a  bounded  response  even  though  the  acoustic  energy  input  is  not  bounded.  This  classic 
nonlinear  behavior  is  attributed  to  the  downward  shift  in  the  column’s  natural  frequency  as  the  deformation 
level  increases.  Under  finite  deformation,  the  natural  frequency  will  be  lower  than  the  value  predicted  by  the 
linear  theory  (Eq.  1);  the  shift  eventually  grows  large  enough  that  the  constant  frequency  acoustic  signal 
begins  to  work  out  of  phase  with  the  column’s  motion.  In  the  case  which  was  investigated,  destructive 
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interference  begins  to  dominate  near  t  =  45  for  the  DRC  injector  design.  Note  that  the  overall  response 
for  the  FRT  is  much  more  chaotic  due  to  the  fact  that  ujn  is  over  seven  times  larger  than  the  acoustic 
frequency  for  this  design.  Rutz5  has  shown  substantial  coupling  in  the  motion  for  frequency  ratios  in  the 
range  0.5  <  ojg/ujn  <  2.  The  fact  that  the  DRC  design  has  a  frequency  ratio  of  1.7  is  the  main  factor 
explaining  its  increased  response  to  the  tangential  mode  frequencies  present  in  the  F-l  engine. 

Figure  5  depicts  the  evolution  of  column  aspect  ratios  to  the  assumed  acoustic  wave  for  both  FRT  and 
DRC  designs.  The  FRT  design  spends  almost  all  the  time  in  the  prolate  mode  (AR  >  1).  This  behavior  is 
due  to  the  fact  that  the  column’s  natural  frequency  is  so  low  (relative  to  the  forcing  frequency)  that  it  is 
responding  primarily  to  the  mean  dynamic  pressure  generated  by  the  acoustic  wave.  The  DRC  design  shows 
a  more  harmonic  response  with  oscillation  between  prolate  and  oblate  shapes.  The  overall  amplitude  of  the 
AR  excursion  is  greater  for  the  DRC  design,  thereby  indicating  more  drastic  changes  in  shape  than  the  FRT 
design. 

This  factor  is  also  evident  in  the  actual  jet  cross-sectional  shapes,  as  revealed  in  Fig.  6.  In  this  figure, 
jet  shapes  are  drawn  to  scale;  the  fact  that  the  DRC  elements  are  smaller  diameter  is  quite  evident.  Here, 
dimensional  times  are  presented  (in  units  of  milliseconds)  as  measured  with  respect  to  the  beginning  of  the 
simulation.  Times  were  selected  to  correspond  with  peaks  in  the  DRC  response  on  the  aspect  ratio  plot.  The 
large  distortions  in  the  DRC  cross-section  are  much  more  pronounced  than  those  of  the  FRT  design.  We 
might  point  out  that  this  conclusion  is  true  in  general;  even  times  at  which  the  FRT  response  is  a  maximum 
reveal  jet  shapes  with  small  distortion  as  compared  to  the  DRC  design. 

Since  the  Weber  number  of  0.1  selected  for  the  previous  simulations  was  arbitrary,  a  series  of  simulations 
were  conducted  to  determine  the  sensitivity  of  both  designs  to  changes  in  this  parameter.  Results  of  these 
simulations  are  presented  in  Fig.  6  in  terms  of  the  actual  peak  speed  in  the  acoustic  wave  ( U  parameter  in 
We).  Here,  simulations  were  run  for  extended  times,  and  the  maximum  aspect  ratio  (in  prolate  form)  was 
determined  for  each  injector.  The  results  indicate  that  the  DRC  design  exhibits  more  deformation  than  the 
FRT  design  at  all  conditions  simulated.  In  fact,  the  deformations  in  the  DRC  design  were  so  large  that  we 
were  unable  to  obtain  stable  solutions  for  peak  acoustic  velocities  exceeding  0.6  m/s  in  this  case.  In  addition, 
the  response  (as  measured  in  terms  of  AR)  is  more  sensitive  to  changes  in  gas  velocity  for  the  DRC  design. 
This  heightened  sensitivity  can  be  viewed  as  a  higher  “gain  factor”  for  pumping  instabilities  in  the  DRC 
design;  i.e.  small  changes  in  velocity  lead  to  rather  substantial  changes  in  jet  distortion. 


Discussion 

There  is  substantial  evidence  that  the  atomization  process  plays  an  important  role  in  many  liquid  rocket 
combustion  instabilities.  Here,  it  is  important  to  note  that  atomization  cannot  in  itself  fully  explain  an 
instability  since  it  is  only  through  resonant  combustion  that  energy  is  made  available  to  enhance  a  small 
perturbation.  However,  basic  atomization  processes  can  be  responsible  for  creation  of  droplet  fields  which 
do  support  resonance.  In  fact,  baffles  are  typically  designed  to  isolate  the  atomization  region  from  the 
combustion  regions11,12. 

The  results  obtained  in  the  previous  section  support  the  conclusion  that  the  frequency  response  of  the 
column  of  fluid  emanating  from  the  orifice  should  be  a  consideration  in  the  design  of  large  liquid  rocket 
engine  injectors.  If  the  column  realizes  substantial  deformation,  aerodynamic  drag  will  be  enhanced  and  the 
impingement  dynamics  will  be  effected.  Early  experimental  work13  indicated  that  liquid  streams  could  in  fact 
respond  to  high-frequency  oscillations.  More  recently,  the  impingement  region14  has  been  shown  to  be  very 
complex  and  highly  dependent  on  alignment  of  the  impact  point.  Depending  on  operating  conditions,  streams 
can  be  reflected,  well  mixed,  or  actually  transmitted  through  the  desired  impact  locale.  Finally,  previous 
experiments  indicate  that  stability  of  impinging  element  injectors  could  be  improved  by  using  contoured 
orifice  inlets11;  the  implication  here  is  that  this  technique  improved  the  accuracy  of  stream  impingement  and 
hence  spray  quality. 

Given  the  sensitivity  of  the  impingement  process  to  jet  alignment,  the  present  scenario  provides  a  rea¬ 
sonable  explanation  for  the  generation  of  instabilities.  Presumably,  droplet  sizes  and  burning  times  are  also 
effected  when  the  impingement  region  atomization  is  varied,  thereby  providing  the  periodic  energy  source 
required  for  instability.  One  could  correlate  the  data  in  Table  1  on  this  basis;  the  frequency  ratio  u)g/ujn  is 
the  key  parameter  in  the  present  scenario.  Figure  8  provides  results  for  this  simple  correlation  for  the  five 
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injector  designs  highlighted  in  Table  1.  Here,  results  indicate  stability  if  the  frequency  ratio  exceeds  a  value 
of  about  3;  beyond  this  value,  the  column  frequency  is  too  low  to  generate  appreciable  response  from  the 
imposed  external  perturbation. 

We  should  note  that  several  attempts  were  made  to  include  data  from  other  testing15^17  and  smaller 
engines12  (such  as  Space  Shuttle  and  Apollo  maneuvering  engines)  with  only  limited  success.  The  bulk  of 
these  tests  utilized  subscale  designs  with  smaller  injector  element  diameters  than  those  which  would  be  used 
in  full  scale  engines.  It  is  possible  that  these  subscale  designs  do  not  provide  sufficiently  intact  jet  boundaries 
upstream  of  the  impingement  point  for  the  present  analysis  to  be  applicable.  In  addition,  it  is  difficult  to 
predict  the  amount  of  acoustic  energy  available  at  the  injector  face  itself;  specular  reflections  within  the 
droplet  field  can  provide  substantial  damping  of  any  presumed  tangential  wave.  Since  the  damping  depends 
on  the  size  and  orientation  of  the  droplet  field,  the  subscale  tests  will  not  replicate  the  full  scale  test  article 
in  this  regard.  Another  factor  involves  the  higher  acoustic  frequencies  associated  with  the  subscale  designs. 
The  length  of  the  acoustic  waves  may  be  so  short  that  a  separation  region  does  not  have  time  to  develop 
on  the  lee-side  of  the  jets,  thereby  limiting  the  deflection  due  to  periodic  drag.  Herein  lies  one  of  the  main 
problems  with  subscale  testing;  while  jet  velocities  can  be  replicated  in  such  tests,  the  atomization  process 
and  resultant  drop  sizes  cannot  be  matched. 

One  should  not  be  so  naive  as  to  expect  a  single  explanation  for  phenomena  as  complex  as  liquid  engine 
combustion  instabilities.  Flame  straining  and  periodic  atomization  theories  can  in  fact  provide  insight 
into  the  overall  process,  but  actual  instabilities  may  occur  as  a  result  of  several  factors.  From  a  practical 
standpoint,  the  F-l  data  suggest  that  larger  orifices  are  more  stable.  This  trend  is  predicted  by  all  the 
present  theories,  including  Hewitt’s  stability  criteria.  However,  the  main  effect  of  increasing  orifice  size  is 
to  move  the  combustion  zone  downstream;  a  design  rule  which  is  known  to  stabilize  many  different  types  of 
injectors.  At  the  present  time,  we  are  forced  to  use  the  limited  theories  as  guidelines  recognizing  that  they 
can  address  very  specific  aspects  of  a  multifaceted  problem. 


Conclusions 

Impinging  element  injectors,  such  as  those  used  on  the  F-l  engine  rely  on  the  injection  of  propellant 
in  well-defined  streams.  These  columns  of  fluid  have  a  natural  frequency  (based  on  linear  analysis)  of 
\/6pg(r/[(p  +  pg)a3].  Nonlinear,  unsteady  boundary  element  simulations  of  the  column  distortion  indicate 
strong  coupling  with  the  acoustic  field  if  the  acoustic  frequency  lies  near  this  column  natural  frequency.  In 
the  case  of  the  F-l  engine,  the  degree  of  instability  (as  measured  by  chamber  pressure  oscillations)  correlates 
well  with  this  frequency  ratio  (u>g/wn).  While  the  actual  F-l  instability  mechanism  could  be  multifaceted, 
this  aerodynamic  interaction  should  be  considered  in  the  design  of  future  large  liquid  engines  utilizing  similar 
impinging  element  injectors. 
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Figure  Captions 

1.  Schematic  of  Impinging-Element  Injector  using  a  Like  Doublet  Configuration.  Tangential  Mode 
Instabilities  Lead  to  a  Periodic  Crosswind  on  Impinging  Fuel  Jets. 

2.  Computational  Domain  and  Boundary  Conditions  for  Simulations  of  Jet  Excitation  from  a  Trans¬ 
verse  Acoustic  Wave 

3.  Grid  Function  Convergence  Test  Results  Showing  Time  History  of  Column  Aspect  Ratio  for 
We  —  0.1,  e  =  0.0064,  u)g/ujn  =  1.7.  Here,  the  Coarse  and  Fine  Grids  have  About  1/2  and  Twice 
the  Number  of  Nodes,  Respectively,  as  in  the  Basis  Grid. 

4.  Motion  of  the  Top  Node  ( x  —  0)  for  Both  DRC  and  FRT  Fuel  Injector  Designs  (e  =  0.0064,  We  - 
0.1). 

5.  History  of  Column  Aspect  Ratio  (Ai?)  for  Both  DRC  and  FRT  Fuel  Injector  Designs  (e  = 
0.0064,  We  =  0.1). 

6.  Fuel-Jet  Cross-Sections  (Plotted  to  Scale)  at  Various  Times  in  Excitation  Process.  Here,  Times 
are  Reported  in  Milliseconds  Relative  to  the  Start  of  the  Simulation. 

7.  Sensitivity  of  Maximum  Jet  Deformation  to  Maximum  Speed  Ug  in  Acoustic  Wave  for  Both  DRC 
and  FRT  Injectors. 

8.  Correlation  of  F-l  Engine  Response  with  Frequency  Ratio,  ojg/ujn. 
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Figure  1:  Schematic  of  Impinging- Element  Injector  using  a  Like  Doublet  Configuration.  Tangential  Mode 
Instabilities  Lead  to  a  Periodic  Crosswind  on  Impinging  Fuel  Jets. 
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Figure  2:  Computational  Domain  and  Boundary  Conditions  for  Simulations  of  Jet  Excitation  from  a  Trans 
verse  Acoustic  Wave 


Dimensionless  Time 


Figure  3:  Grid  Function  Convergence  Test  Results  Showing  Time  History  of  Column  Aspect  Ratio  for 
We  =  0.1,  e  =  0.0064,  0Jg/uJn  =  1.7.  Here,  the  Coarse  and  Fine  Grids  have  About  1/2  and  Twice  the 
Number  of  Nodes,  Respectively,  as  in  the  Basis  Grid. 
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Figure  6:  Fuel-Jet  Cross-Sections  (Plotted  to  Scale)  at  Various  Times  in  Excitation  Process.  Here,  Times 
are  Reported  in  Milliseconds  Relative  to  the  Start  of  the  Simulation. 
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7  Appendix  C  -  Hybrid  Rocket  Injection/ Combustion  Tests 


Wernimont,  E.J.,  and  Heister,  S.D.,  “Combustion  Experiments  in  a  Hydrogen  Perox¬ 
ide/Polyethylene  Hybrid  Rocket  with  Catalytic  Ignition”,  In  Review,  Journal  of  Propulsion 
and  Power ,  1998. 
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COMBUSTION  EXPERIMENTS  IN  A  HYDROGEN 
PEROXIDE/POLYETHYLENE  HYBRID  ROCKET 
WITH  CATALYTIC  IGNITION 


E.  J.  Wernimont  *and  S.  D.  Heister  t 
Purdue  University,  West  Lafayette,  Indiana  47907 

January  28,  1999 


Abstract 

Results  from  100  tests  of  labscale  hybrid  rocket  motors  using  hydrogen  peroxide  and  polyethy¬ 
lene  are  presented  herein.  The  bulk  of  the  tests  utilized  85%  peroxide  with  low  density  polyethy¬ 
lene.  A  new  consumable  catalytic  ignition  device  has  been  utilized  to  provide  rapid,  reliable 
ignition  using  stabilized  peroxide.  Regression  measurements  indicate  that  at  low  chamber  pres¬ 
sures  (100  psi)  a  classic  diffusion-dominated  behavior  is  noted  with  massflux  exponents  very 
near  the  theoretical  value  of  0.8.  However,  at  higher  chamber  pressures  tested  (200  and  400 
psi),  radiative-dominated  behavior  is  noted  for  average  massfluxes  varying  between  0.1  and  0.3 
lbm/(in2  s).  Through  the  optimization  of  aft  mixing  length,  combustion  efficiencies  in  excess  of 
95%  were  obtained  in  these  tests.  No  significant  nonacoustic  or  acoustic  instabilities  were  noted 
in  these  tests;  chamber  pressure  fluctuations  were  less  than  3.5%  zero-to-peak  of  mean. 


Nomenclature 

At  Throat  area  (in2) 

C*  Characteristic  exhaust  velocity  (ft/s) 

dt  Differential  time  (sec) 

Dp  Port  Diameter  (in) 

G  Total  mass  flux  (lbm/(in2s)) 

gc  Gravitational  constant  (^yfr) 

L j  Fuel  grain  length  (in) 

L *  Characteristic  chamber  dimension  (in) 

Mf  Fuel  mass  (lbm) 

Mox  Oxidizer  mass  (lbm) 

Mccb  CCB  mass  (lbm) 

Minert  Mass  of  expended  inerts  (lbm) 
mox  Oxidizer  mass  flow  rate  (lbm/s) 

n  Time  level 

OF  Mass  oxidizer/fuel  mixture  ratio 

Pc  Head-end  chamber  pressure  (psia) 

r  Fuel  regression  rate  (in/s) 

t  Time  (sec) 

*  Chief  Propulsion  Engineer,  Beal  Aerospace,  Dallas,  TX 

*  Associate  Professor,  School  of  Aeronautics  and  Astronautics.  Associate  Fellow,  AIAA 
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pj  Fuel  density  (lbm /in3) 


Super  /  Subscript 

()  Burning-time  averaged  quantity 

f  Final 

i  Initial 

Acronyms 

CCB  Consumable  Catalytic  Bed 

HDPE  High  Density  Polyethylene 

HP  Hydrogen  Peroxide 

LDPE  Low  Density  Polyethylene 
PMMA  Polymethyl  Methacrylate 
UHMW  Ultra  High  Molecular  Weight 


Introduction 

In  the  past  few  years,  interest  in  hybrid  rockets  has  increased  due  to  the  potential  for  these  devices  to  reduce 
costs  and  enhance  safety  in  aerospace  propulsion  devices.  A  variety  of  applications,  including  launch  vehicle 
boosters,  upper  stage  and  tactical  systems  have  been  identified  as  areas  in  which  hybrid  propulsion  concepts 
are  of  interest. 

We  can  trace  the  use  of  the  Hydrogen  Peroxide  (HP)/  Polyethylene  (PE)  propellant  combination  to  the 
mid  1950’s1”3.  While  the  early  tests  of  Moore  and  Berman1  were  quite  successful,  interest  waned  (most 
probably  due  to  the  search  for  higher  energy  propellants  during  this  era)  and  essentially  no  published  work 
exists  for  a  four  decade  period  beginning  in  the  mid  1950’s.  Current  requirements  for  lower  cost,  non-toxic 
propulsion  systems  have  motivated  a  renewed  interest  in  this  storable  propellant  combination.  Recent  efforts 
have  been  undertaken  in  our  group4”9,  at  the  University  of  Surrey  in  the  United  Kingdom10  and  the  U.S. 
Air  Force  Academy11. 

Recent  system  studies4,5,12  point  to  potential  benefits  of  HP-oxidized  systems  which  include  its  high 
density,  ease  of  handling,  non-toxicity,  and  monopropellant  characteristics.  For  example,  both  turbopump 
and  pressurization  systems  can  utilize  decomposition  energy  and  biproducts  to  effectively  simplify  engine 
power  cycles  and  tank  pressurization  systems.  In  addition,  the  fact  that  HP  hybrid  systems  tend  to  optimize 
at  high  mixture  ratios  provides  a  benefit  in  reducing  the  size  of  the  expensive,  high-pressure  combustion 
chamber  which  contains  the  fuel  grain.  This  benefit  also  decreases  the  sensitivity  of  these  propellants  to  fuel 
slivering  since  the  fuel  provides  a  smaller  fraction  of  the  total  propellant  mass. 

Polyethylene  (PE)  also  represents  a  unique  fuel  choice  in  view  of  the  fact  that  much  of  the  present  work 
focuses  on  the  use  of  hydroxyl-terminated  polybutadiene  (HTPB)  as  fuel.  While  HTPB  systems  enjoy  a 
slightly  higher  regression  rate,  which  tends  to  reduce  fuel  grain  complexity  somewhat,  this  fuel  is  a  thermoset 
compound  which  can  only  be  produced  using  batch  processing  which  necessitates  substantial  tooling  capital 
investment,  tooling  assembly  and  disassembly.  We  suspect  that  the  main  impetus  for  the  use  of  this  chemical 
stems  from  the  fact  that  most  solid  rocket  manufacturers  (many  of  whom  also  participate  in  hybrid  rocket 
programs)  commonly  use  this  chemical  in  solid  rocket  propellants. 

In  contrast  to  HTPB,  PE  is  a  thermoplastic  which  is  commonly  produced  via  extrusion  from  a  die  in  a 
continuous  process.  Hence,  PE  grains  could  be  produced  using  this  approach  by  simply  cooling  the  extruded 
product  and  cutting  it  to  the  desired  length.  Thermoplastics  also  eliminate  waste  since  the  product  can  be 
remelted  if  a  part  is  made  incorrectly.  In  addition,  PE  has  a  lower  cost  than  HTPB  and  is  much  easier  to 
machine  even  though  the  thermochemical  combustion  performance  of  the  two  materials  is  virtually  identical. 
For  these  reasons,  this  material  represents  an  attractive  alternative  for  many  missions. 

For  these  reasons,  an  experimental  program  was  initiated  to  quantify  the  combustion  behavior  of  the 
HP/PE  hybrid  propellant  combination.  A  unique  consumable  catalytic  ignition  system  was  used  to  provide 
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the  initiation  of  combustion  in  these  studies.  Factors  considered  in  the  tests  described  in  this  paper  include: 
massflux,  oxidizer/fuel  (OF)  mixture  ratio,  chamber  pressure,  fuel  grain  length  ( L *),  and  PE  formulation. 
The  facility  is  briefly  described  in  the  following  section,  followed  by  a  description  of  the  ignition  system,  and 
experimental  results. 

Test  Apparatus/Methodology 

This  fluid  system  was  designed  to  provide  reliable  combustion  measurements  with  minimal  complexity  and 
redundant  safety  features.  To  this  end,  a  simple  cartridge-loaded  combustion  chamber  is  utilized  with  a 
simple  “nozzle”  designed  to  provide  the  required  throat  area  without  an  exit  cone  (simple  throat  plug). 
Figure  1  highlights  the  2  inch  outside  diameter  combustion  chamber  design;  the  catalytic  ignition  device  is 
described  in  detail  in  the  following  section.  The  insulating  materials  used  for  these  tests  are  paper  phenolic. 
The  nozzle  material  is  a  silica  phenolic  that  provided  a  low  erosion  rate  of  less  than  2.0  mils/sec  for  the 
configurations  fired. 

A  schematic  of  plumbing  and  tankage  associated  with  the  test  apparatus  is  shown  in  Fig.  2.  This  appa¬ 
ratus  provides  for  safe  firing  of  a  hydrogen  peroxide  hybrid  rocket  motor  by  remote  operation  from  a  concrete 
enclosed  control  room.  The  entire  fluid  system  utilizes  materials  compatible  with  high-concentration  hydro¬ 
gen  peroxide  such  as  300-series  stainless  steels,  glass  and  teflon.  The  high  pressure  HP  tank  is  pressurized 
using  a  Nitrogen  “K-bottle”  with  manual  isolation  valve,  MV1,  depicted  in  Figure  2.  A  second  K-bottle 
(manual  isolation  valve  MV4)  is  used  to  provide  high  pressure  gas  for  remote  actuation  of  stainless  steel 
pneumatic  valves  denoted  PV1  and  PV2  in  Fig.  2.  During  a  test,  the  system  is  operated  by  simply  opening 
the  main  oxidizer  valve  (PV1)  until  all  of  the  HP  loaded  in  the  oxidizer  tank  is  consumed.  Gaseous  nitrogen 
is  then  driven  into  the  combustion  chamber  causing  quenching  of  the  remaining  fuel. 

As  can  be  seen  from  Fig.  2,  this  system  incorporates  many  safety  features.  Actuation  of  pneumatic  valve 
PV2  provides  for  dumping  and  dilution  of  the  HP  (PV2)  in  the  event  of  an  emergency.  During  the  entire  test 
program,  there  were  no  occasions  in  which  dumping  of  the  HP  was  required.  Other  safety  features  include: 
a  normally-open  venting  solenoid  valve  (SV2),  a  remotely-actuated  manual  vent  valve  (MV5),  and  a  relief 
valve  (RV1)  on  the  oxidizer  tank.  These  features  were  installed  to  allow  venting  of  the  oxidizer  system  in 
the  event  that  undesired  decomposition  of  the  loaded  HP  occurs. 

Instrumentation  for  these  tests  consist  of  ullage  and  chamber  pressure,  oxidizer  turbine  flow  meter 
information  as  well  as  axial  thrust  measurement.  This  set  of  information  (plus  pre/post-test  fuel  grain 
measurements)  is  sufficient  to  determine  motor  operating  parameters  as  well  as  information  on  fuel  regression 
rate.  Digital  data  acquisition  is  achieved  through  the  use  of  a  Pentium-based  PC  with  a  Keithley/Metrabyte 
analog/digital  converter.  This  data  acquisition  system  is  used  to  real  time  display  the  conditions  of  oxidizer 
tank  temperature  (Tl)  and  pressure  (P3)  during  HP  loading.  This  provides  a  means  of  determining  if  an 
emergency  dumping  procedure  must  be  conducted  if  the  loaded  HP  is  unstable  in  the  tank. 

During  a  firing  the  data  acquisition  system  is  used  to  sample  the  motor  operation  parameters  at  250 
Hz  per  instrument.  In  each  test,  oxidizer  flow  rate  is  held  approximately  constant  using  the  self-limiting 
combustion  behavior  of  hybrid  rocket  motors,  i.e.,  the  rate  of  oxidizer  injected  into  the  combustion  chamber 
is  a  direct  function  of  the  ullage-to-ch amber  pressure  differential  and  the  chamber  pressure  is  also  dependent 
on  the  oxidizer  flowrate  (for  given  nozzle  and  fuel  port  diameters).  Assuming  no  throat  erosion,  both  chamber 
pressure  and  oxidizer  flow  rate  may  be  held  constant  by  maintaining  a  constant  ullage  pressure.  Constant 
ullage  pressure  is  achieved  by  regulated  nitrogen  pressurant  gas.  This  technique  is  used  for  all  the  tests 
presented  in  this  paper. 

The  test  program  was  formulated  to  maximize  the  amount  of  information  available  to  actually  design 
a  HP/PE  hybrid.  For  this  reason,  many  different  parameters  were  studied,  none  of  which  were  studied  in 
tremendous  detail.  Most  parameters  were  investigated  across  a  range  of  total  mass  port  flux  levels  allowing 
determination  of  fuel  regression  rate  influence.  The  entire  test  program  involved  firing  eleven  different  series 
of  motors,  with  each  series  dedicated  to  examination  of  a  specific  parameter.  A  total  of  100  firings  were 
conducted  in  association  with  these  test  series  which  are  summarized  in  Table  1.  Acronyms  in  Table  1  are 
defined  in  the  nomenclature.  The  following  sections  will  discuss  results  obtained  from  the  bulk  of  these  tests. 
Testing  accomplished  using  a  radial-flow  geometry  (test  series  C  in  Table  1)  are  reported  in  Ref.  13. 
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Table  1:  Test  Series  Summary 


Series 

Purpose 

Number 
of  Motors 

Test  Parameters 

A 

System  Verification 

12 

Pc  ~  100  psia 

G  0.25  to  0.35  lbm / (in2  —  s) 

85%  HP  Interox,  HDPE 

B 

Broader  Flux 

12 

G  0.15  to  0.35  lbm/(in2  —  s) 
Three  at  80%  HP  Interox 

C 

Radial  Flow  Motor 

20 

Proof  of  Concept 

8  Ignition  Tests 

D 

Polyethylene  Types 

11 

LDPE,  UHMW 

E 

Higher  Mass  Flux,  Pc 

8 

Flight  Weight,  LDPE 

Proof  of  Combustion 

F 

New  HP  Vendor 

5 

Air  Liquide  85%  HP 

G 

Increase  Pc 

5 

Pc  ~  200  psia 

H 

Lengthen  Aft  Combustion 

4 

Increase  by  2,  4  Inches 

I 

Action  Time  Study 

4 

Two  at  6,  9  Seconds  Each 

J 

Increase  Pc 

5 

Pc  ~  400  psia 

M 

Increase  Mass  Flux 

10 

G  0.35  to  1.0  lbm/(in2  -  s) 

Y 

Ignition  Behavior 

4 

Visual  Observations,  PMMA  Fuel 

Consumable  Ignition  Device 

Probably  the  greatest  challenge  in  creating  a  workable  hybrid  propellant  combination  lies  in  the  development 
of  an  ignition  concept  which  provides  a  rapid  and  reproducible  rise  in  chamber  pressure  and  thrust.  In  the 
past,  secondary  injection  of  pyrofuoric  fluids,  electric  ignition  sources,  torches,  and  catalytic  ignition  systems 
have  been  used  in  hybrid  rockets.  With  the  exception  of  the  catalytic  system,  all  these  concepts  require 
additional  hardware  and/or  fluids  to  support  ignition  of  the  motor.  For  this  reason,  the  catalytic  concept 
was  pursued  through  the  use  of  a  consumable  catalytic  bed  (CCB).  While  catalytic  devices  such  as  silver 
screens  and  other  materials  treated  with  catalytic  material  have  been  used  to  decompose  peroxide  in  the 
past14,  the  present  concept15  provides  several  advantages  over  these  techniques: 

•  The  CCB  concept  supports  operation  with  stabilized  HP  (successful  tests  with  stabilizer  levels  as  high 
as  50  ppm).  Most  of  the  previous  systems  could  not  operate  with  stabilizers  (primarily  stannate  and 
phosphate  compounds)  in  the  fluid  because  these  materials  would  reduce  the  catalytic  activity  of  the 
bed  material.  Use  of  stabilized  fluid  throughout  this  test  program  has  been  viewed  as  a  substantial 
safety  enhancement. 

•  The  CCB  provides  an  energetic  ignition  system  with  no  inert  parts,  thereby  maximizing  the  propellant 
mass  fraction  of  the  propulsion  system. 

•  Since  the  device  requires  no  inert  hardware  (or  other  fluids/gases)  it  also  represents  a  minimum  cost 
ignition  system. 

The  CCB  is  inserted  into  a  pocket  which  was  machined  into  the  forward  end  of  the  fuel  grain  as  noted 
schematically  in  Fig.  L  No  special  ignition  sequence  was  utilized  in  any  of  the  testing,  the  engine  was 
literally  “slam  started”  by  opening  valve  PV2  (see  Fig.  2)  which  provides  maximum  HP  flow.  As  HP  passes 
through  passages  in  the  CCB  it  is  decomposed.  If  the  HP  is  at  sufficient  concentration,  the  decomposition 
products  are  at  a  temperature  such  that  autoignition  of  the  PE  occurs.  As  the  HP  flow  continues,  the  CCB 
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is  consumed.  If  the  CCB  is  sized  properly,  there  will  be  enough  energy  in  the  combusting  fuel  grain  to 
support  thermal  decomposition  of  the  HP  injected  after  CCB  consumption;  i.e.  the  device  operates  for  only 
a  small  fraction  of  a  firing. 

Figure  3  highlights  a  group  of  typical  ignition  transients  obtained  with  the  use  of  85%  HP  in  the  design. 
This  collection  shows  that  the  CCB  is  repeatable  and  produces  rapid  ignition  transients.  The  peaks  at 
approximately  30  msec  represent  a  pressure  overshoot  due  to  increased  oxidizer  mass  flow  when  the  oxidizer 
valve  is  first  opened.  Because  there  is  no  cavitating  venturi  in  the  fluid  system,  the  oxidizer  mass  flow  is 
governed  strictly  by  pressure  differential.  Consequently,  initial  fluid  flow  is  greatest  for  the  first  20  msec 
before  the  chamber  pressure  has  a  chance  to  climb  to  its  steady  state  value.  As  can  be  seen  from  the  figure, 
steady-state  combustion  has  been  achieved  in  50  msec. 

To  aid  in  the  understanding  of  the  rate  of  consumption  of  the  device,  a  series  of  four  motors  (Y-series 
described  in  following  section)  were  fired  using  transparent  acrylic,  (PMMA)  fuel  grains.  A  typical  chamber 
pressure  history  from  one  of  these  firings  (test  Yl)  is  shown  in  Fig.  4.  The  time  window  in  the  figure 
captures  the  interval  from  ignition  to  the  beginning  of  tailoff  such  that  the  ignition  event  occurs  at  roughly 
3.3  sec  on  this  time  scale.  There  is  a  noticeable  decrease  in  chamber  pressure  at  roughly  5.5  sec  (2.1  sec 
after  ignition)  which  is  observed  in  most  firings  in  our  test  program.  This  behavior  was  attributed  to  the 
complete  combustion  of  the  CCB.  A  similar  behavior  was  noted  on  other  PMMA  firings  and  rapid  decreases 
in  chamber  pressure  were  well  correlated  with  consumption  of  the  CCB. 

To  investigate  this  assertion,  still  photographs  were  taken  during  this  firing.  The  photos  were  taken  with 
the  auto-iris  set  on  the  camera  so  intensity  levels  are  not  necessarily  representative  of  flame  temperature.  A 
sequence  of  four  photos  is  summarized  in  Fig.  5.  In  this  figure,  the  CCB  is  on  the  right  and  flow  is  from 
right  to  left.  At  ignition  (t=3.3  sec),  the  fuel  port  is  already  luminous  indicating  reaction  within  the  CCB. 
The  second  image,  taken  just  prior  to  the  chamber  pressure  decrease,  shows  complete  consumption  of  a  small 
portion  of  the  CCB  at  one  circumferential  location.  In  the  next  frame  available  from  the  camera  (t=5.5  sec) 
the  entire  aft  portion  of  the  CCB  has  been  consumed  or  expelled  from  the  motor.  This  event  correlates  well 
with  the  decrease  in  chamber  pressure  noted  in  Fig.  4.  The  final  image,  taken  just  prior  to  total  oxidizer 
consumption  (t=11.5  sec)  shows  a  small  dark  region  at  the  head-end  of  the  grain.  Since  postfire  inspection 
revealed  that  the  entire  CCB  had  been  consumed  (or  expelled),  this  dark  region  corresponds  to  the  zone  in 
which  HP  decomposition  occurs.  Postfire  fuel  samples  are  in  agreement  with  this  theory;  negligible  regression 
is  observed  in  this  region. 

One  can  also  note  an  asymmetric  fuel  regression  pattern  on  the  last  image  (t=11.5  sec);  maximum  fuel 
regression  occurs  at  the  top  of  the  CCB  pocket  in  a  region  which  was  the  last  area  to  be  exposed  a  result  of 
local  CCB  consumption.  This  asymmetric  behavior  was  noted  in  all  four  PMMA  firings  and  was  attributed  to 
nonuniformities  in  the  spray  pattern  emanating  from  the  injector16.  The  higher  flow  which  caused  local  CCB 
consumption  at  t=5.45  sec.  could  provide  local  quenching  (or  decreased  combustion),  thereby  explaining 
the  observed  behavior.  While  some  asymmetries  were  noted  in  firings  using  PE  grains,  they  were  generally 
restricted  to  the  region  which  housed  the  CCB. 

Test  Results 

As  indicated  in  Table  1,  a  broad  range  of  tests  were  conducted  over  the  twelve  test  series.  Unfortunately, 
we  were  unable  to  maintain  a  consistent  HP  vendor  over  the  three  year  study.  While  some  differences  in 
fluid  were  noted16,  they  were  generally  minor.  In  addition,  we  experimented  with  several  different  throat 
materials  before  arriving  at  a  suitable  option,  which  was  used  throughout  test  series  D-Y.  Early  efforts  were 
aimed  at  optimizing  fuel  grain  length  so  as  to  obtain  OF  ratios  near  the  7.5  value  which  tends  to  maximize 
specific  impulse  for  these  propellants.  Many  of  the  later  tests  actually  achieved  fairly  low  mixture  ratios  (in 
the  5-6  range)  due  to  the  fact  that  fuel  regression  rates  exceeded  our  estimates.  A  complete  tabulated  list  of 
all  test  conditions  can  be  found  in  Ref.  16;  we  will  not  include  all  these  data  here  in  the  interest  of  brevity. 

Substantial  efforts  were  expended  in  quantifying  combustion  behavior  for  various  PE  formulations  in¬ 
cluding  Low  Density  PE  (LDPE),  High  Density  PE  (HDPE),  and  ultra-high  molecular  weight  (UHMW) 
materials16.  Discussion  of  this  area  is  planned  for  a  future  work;  the  bulk  of  the  efforts  to  be  described 
herein  relate  to  the  effects  of  chamber  pressure,  fuel  grain  internal  diameter  and  length,  aft  combustion 
length  and  OF  ratio  on  combustion  performance. 

Because  the  CCB  device  contributed  a  non-negligible  amount  of  mass  and  energy  during  our  test  burns, 
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which  typically  varied  between  5  and  20  seconds,  we  were  forced  to  improve  upon  the  data  reduction  method¬ 
ology  which  had  typically  been  used  by  hybrid  rocket  experimentalists.  Average  regression  rates  and  port 
total  flux  levels  from  a  given  test  firing  are  obtained  from  an  integral  reconstruction  of  the  entire  cham¬ 
ber  pressure  history  with  the  CCB  and  inerts  flows  taken  into  account.  The  approach  begins  with  the 
determination  of  average  characteristic  velocity  (C*)  and  oxidizer/fuel  mass  ratio  (OF)  for  a  given  test: 

=  /;/  Pc(t)At(t)gedt 

Mox  +  Mj  -f-  Mcch  +  Minert 


and 


OF  = 


MqX 

Mfuei  +  Mccb  +  Minert 


(2) 


The  expended  masses  (see  Nomenclature  for  definitions)  can  all  be  measured  directly  with  the  use  of  pre-  and 
postfire  measurements.  Instantaneous  values  (time  level  n)  of  total  mass  flux  (from  its  definition)  and  fuel 
regression  rate  (assuming  a  steady-state  balance  of  incoming  and  exiting  mass  flows)  may  then  be  calculated 
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Using  this  process,  the  entire  regression  rate  and  mass  flux  histories  may  be  reconstructed  from  the  measured 
time-dependent  data.  Action  time  averages  are  computed  via  direct  integration: 


ftf 

r  =  J  rn(t)dt/  J  dt 

(5) 

ftf  ftf 

(6) 

G=  Gn(t)dt/  dt 

Jtt  Jti 

Error  analyses  indicate  that  for  our  apparatus,  this  approach  is  superior  to  more  classical  “end-point”  based 
techniques17  which  compute  average  regression  rate  and  port  total  mass  flux  using  just  pre-  and  postfire 
port  diameters.  For  the  tests  presented  herein,  the  bias  errors  in  regression  rate  and  flux  are  estimated  to  be 
1.9  and  1.1%  respectively.  A  complete  description  of  the  methodology,  which  will  be  used  throughout  this 
work,  is  provided  in  Ref.  17. 

Reliable  ignition  and  combustion  was  demonstrated  over  a  range  of  initial  oxidizer  fluxes  0.1  <  Gox  < 
1.2  lbm/(in2  s)  and  chamber  pressures  100  <  Pc  <  400  psia  during  the  testing  with  85%  HP.  A  few  tests 
were  conducted  with  80%  HP,  but  reliable  ignition  and  combustion  could  not  be  achieved  using  the  injector 
and  CCB  design  implemented  in  these  tests.  We  believe  that  one  could  design  a  device  to  operate  at 
these  lower  concentrations  with  an  improved  injector  (with  smaller  droplet  sizes)  and  CCB  designs.  Our 
injector  produced  droplets  which  were  quite  large,  a  350  micron  volumetric  mean  diameter  is  reported  by 
the  manufacturer.  Limited  information16  suggests  that  ignition  and  combustion  at  concentrations  lower 
than  80%  would  be  very  difficult  for  PE;  at  these  concentrations  most  of  the  decomposition  energy  goes  into 
vaporizing  the  water  in  the  aqueous  HP.  Here  we  should  note  that  concentrations  below  67%  have  insufficient 
decomposition  energy  to  vaporize  water  within  the  mixture. 

Typical  chamber  pressure  histories  obtained  during  the  test  program  are  shown  in  Fig.  6.  Overall, 
the  combustion  obtained  was  very  smooth  as  evidenced  by  the  low  amount  of  noise  in  the  pressure  signals. 
Typical  unsteadiness  in  the  pressure  signals  was  of  the  order  of  1-2%  (zero-to-peak)  of  the  mean  pressure  in 
this  test  program.  Sharp  tailoffs  were  always  observed  in  the  testing;  action  times  were  determined  using  a 
bisector  technique17  commonly  used  in  solid  rocket  data  reduction.  The  spike  in  the  motor  F4  and  H3  traces 
in  the  interval  1  <  t  <  2  seconds  is  attributed  to  ejection  of  small  portions  of  the  CCB  through  the  nozzle. 
In  the  Ml  pressure  trace,  a  slight  increase  in  Pc  is  observed  after  the  main  ignition  event.  This  behavior 
may  be  attributed  to  increased  oxidizer  flow  rate  due  to  a  change  in  injector  discharge  coefficient16. 
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Table  2:  Summary  of  LDPE/85%  HP  Fuel  Regression  Flux  Correlations 


Correlation 

Pc 

G 

0.040G0'78  ips 
0.035G0-52  ips 
0.041G0'49  ips 

100  psia 
200  psia 
400  psia 

0.1  -  0.3  lbm/(in2  —  s) 
0.1 -0.3  lbm / (in2  —  s) 
0.2 -0.7  lbm/(in2  —  s) 

Regression  Rate  Behavior 

While  there  were  test  series  dedicated  to  assessing  the  influence  of  PE  type  on  fuel  regression  and  combustion, 
the  bulk  of  the  measurements  were  obtained  using  Low  Density  PE  (LDPE)  fuel.  A  compilation  of  these 
measurements  is  provided  in  Fig.  7  which  highlights  dependence  of  regression  rate  on  both  flux  level  and 
chamber  pressure.  While  the  low  pressure  (100  psi)  data  appear  to  behave  a  classical  regression  law  (r  oc  Gn) 
consistent  with  diffusion-dominated  behavior,  the  higher  pressure  results  show  a  distinct  insensitivity  to  port 
total  mass  flux  (G)  in  the  range  0.1  <  G  <  0.3  lbm/(in2  s).  In  this  lower  flux  region,  regression  rates 
at  the  higher  pressures  tested  (200  and  400  psia)  are  as  much  as  75%  greater  than  those  at  low  pressures. 
The  data  are  consistent  with  a  radiation-based  regression  law  which  has  been  theorized,  but  infrequently 
observed  for  low  massflux  conditions. 

In  fact,  this  behavior  is  ideal  since  regression  rates  would  no  longer  be  influenced  by  changes  in  port 
geometry  (either  shape  or  size).  For  the  booster  application,  design  studies5  indicate  an  optimal  massflux 
level  of  about  0.4  lbm/(in2  s)  assuming  a  classical  regression  law,  r  oc  G°  8.  Presumably,  designers  could 
make  use  of  this  desirable  behavior  (flux  insensitivity)  for  other  applications  as  well. 

To  compare  the  data  from  this  test  program  with  that  of  previous  researchers,  we  performed  correlations 
assuming  a  classical,  massflux-dominated  regression  behavior.  While  this  approach  is  not  warranted  for  the 
higher  pressure  results,  it  does  permit  gross  comparisons  with  results  of  other  researchers  who  made  similar 
assumptions.  The  resulting  correlations  are  presented  in  Table  2.  Note  that  for  the  low  pressure  data  the 
exponent  of  0.78  is  in  close  agreement  with  theory  assuming  a  diffusion-dominated  behavior.  The  exponent 
is  reduced  at  the  higher  pressures  due  to  the  radiative-dominated  behavior. 

The  correlations  from  Table  2  are  compared  with  those  obtained  by  other  researchers  in  Fig.  8.  Previous 
researchers  using  hydrogen  peroxide1-3,18,19  had  also  noted  a 

Figure  9  shows  combustion  efficiency  versus  the  length  of  the  aft  combustion  chamber.  As  can  be  seen 
from  the  figure,  there  appears  to  be  significant  efficiency  gains  from  even  minor  length  increases.  There  are 
only  two  motors  fired  with  a  6  inch  length  since  it  is  thought  sufficient  gains  are  made  firing  with  a  4  inch 
long  aft  combustion  chamber.  Note  that  efficiency  values  of  above  95%  are  routinely  attained  in  this  device 
with  minor  changes  in  combustor  geometry. 

Figure  10  shows  the  combustion  efficiency  versus  the  characteristic  chamber  length,  L*,  for  the  motors 
tested.  For  our  purposes  L*  is  calculated  as  the  combustor  volume  aft  of  the  fuel  grain  and  forward  of  the 
throat  divided  by  the  throat  area.  This  parameter  is  often  used  in  the  liquids  industry  to  design  a  motor 
for  acceptable  (>95%)  combustion  efficiency.  This  parameter  is  specific  to  propellant  combination,  but  for 
comparison  to  liquid  oxygen  and  ethyl  alcohol  L*  values  are  between  40  and  120  inches21.  This  figure  shows 
that  efficiency  gains  are  derived  from  increased  L*  as  is  generally  the  case.  This  figure  suggests  that  L* 
values  as  low  as  40  inches  may  provide  sufficient  aft  combustion  volume  for  efficiencies  above  90%. 

Figure  11  shows  the  effect  of  measured  average  chamber  pressure  on  combustion  efficiency.  This  figure 
also  tends  to  suggest  that  combustion  efficiency  gains  are  obtained  as  a  result  of  increases  in  chamber 
pressure  and  that  values  of  at  least  200  psi  are  required  to  obtain  acceptable  efficiencies  (at  least  for  the  aft 
combustion  lengths  tested).  One  could  speculate  that  the  100  psi  firings  may  have  required  a  larger  mixing 
region  to  obtain  a  higher  efficiency;  this  theory  was  not  explored  in  our  testing.  Since  most  applications 
demand  chamber  pressures  in  excess  of  100  psi,  this  issue  should  not  adversely  effect  a  given  design. 
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Combustion  Stability  Behavior 

As  noted  in  the  discussion  of  Fig.  6,  there  were  no  notable  combustion  instabilities  in  this  test  program. 
Recent  theoretical  efforts22  suggest  that  thermal  lags  in  the  fuel  can  account  for  the  low-frequency  (1-100  Hz) 
oscillations  which  have  been  observed  in  numerous  tests  which  have  used  either  liquid  or  gaseous  oxygen 
as  the  oxidizer.  In  contrast  to  liquid  oxygen  systems  which  tend  to  optimize  at  OF  values  near  2.5,  the 
HP/PE  system  optimizes  at  much  higher  OF  ratio  (typically  around  7.5).  In  this  case,  the  mass  and  energy 
addition  from  the  fuel  is  much  less  than  in  oxygen  oxidized  systems,  thereby  reducing  the  overall  amplitude 
achievable.  This  factor  may  explain  our  distinct  lack  of  instabilities.  However,  we  should  note  that  all  testing 
was  performed  on  a  fairly  small  scale  apparatus  (2  inch  diameter  fuel  grain)  and  additional  larger  scale  work 
will  be  required  to  confirm  this  suspicion. 

Overall,  the  level  of  combustion  oscillations  was  observed  to  be  below  3.5%  zero-to-peak  of  mean  through¬ 
out  the  testing  program.  Figure  12  shows  a  typical  waterfall  plot  from  motor  F4.  A  slight  tendency  for 
non-acoustic  combustion  oscillations  up  to  a  frequency  of  approximately  70  Hz  is  noted.  Comparison  with 
behavior  in  other  motors  suggests  that  the  higher  frequency  (acoustic)  portion  may  account  for  roughly  1.5% 
zero-to-peak  chamber  pressure  response.  Although  the  motor  F4  waterfall  plot  shows  no  preferred  waveform 
for  this  propellant  combination,  later  motors  do  have  preferred  oscillatory  waveforms.  As  an  example  Fig¬ 
ure  13  shows  the  waterfall  plot  for  motor  G3  which  has  a  preferred  waveform  at  roughly  20  to  35  Hz  with  an 
amplitude  of  around  0.75  psid  zero-to-peak.  Three  other  motors  (G4, 14  and  M3)  exhibited  similar  responses 
to  that  of  motor  G3.  Other  motors  (J3,  G5,  and  J5)  exhibited  a  preferred  waveform  at  roughly  65  to  80  Hz 
with  an  amplitude  of  around  0.5  psid  zero-to-peak. 

Conclusions 

This  paper  summarizes  combustion  measurements  from  testing  of  the  hydrogen  peroxide  (HP) /polyethylene 
(PE)  hybrid  rocket  propellant  combination.  The  bulk  of  the  tests  utilized  85%  HP  with  low  density  PE 
(LDPE),  but  some  data  are  reported  for  other  types  of  PE.  A  new  consumable  catalytic  bed  (CCB)  design 
has  been  used  to  provide  rapid,  reproducible  ignition  using  stabilized  HP.  Once  the  CCB  is  consumed,  the 
HP  undergoes  thermal  decomposition  as  a  result  of  exposure  to  combustion  gases  emanating  from  the  ignited 
fuel  grain.  The  device  provides  a  simple,  low  cost  (and  weight)  alternative  for  ignition  of  hybrid  rockets. 

Measured  regression  rates  indicate  a  classic  diffusion-dominated  behavior  at  chamber  pressures  of  100  psi, 
with  flux  exponents  very  near  the  theoretical  value  of  0.8.  However,  at  higher  pressures,  radiation-dominated 
behavior  in  which  regression  appears  to  be  insensitive  to  changes  in  flux,  is  noted  for  massflux  levels  between 
0.1  and  0.3  lbm/(in2  s).  This  behavior  can  lead  to  pressure-related  amplifications  of  regression  rate  as  high 
as  75%  when  compared  to  the  low  pressure  (100  psi)  behavior.  Results  are  shown  to  be  comparable  to  those 
obtained  from  other  researchers  using  HP  as  oxidizer  with  various  other  fuels. 

High  combustion  efficiencies  (>  95%)  were  obtained  at  the  higher  chamber  pressure  (200  and  400  psi) 
conditions  by  using  aft  mixing  lengths  equivalent  to  about  two  fuel  grain  diameters  (4  inches) .  These  lengths 
are  consistent  with  L*  values  (based  on  the  chamber  volume  aft  of  the  fuel  grain)  of  about  60  inches. 
Smooth  combustion  was  observed  in  all  testing,  with  typical  chamber  pressure  fluctuations  in  the  range  of 
1-2%  (zero-to-peak)  of  the  mean.  We  speculate  that  the  high  mixture  ratio  (5-8)  operation  of  this  propellant 
combination  plays  a  role  in  reducing  the  amount  of  energy  available  to  drive  nonacoustic  instabilities.  Some 
minor  activity  is  noted  in  the  range  of  20-80  Hz  for  some  of  the  100  motors  which  were  tested  in  these  efforts. 
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Figure  Captions 

1.  Cross  Sectional  View  of  the  Combustion  Chamber 

2.  Test  Facility  Schematic 

3.  Typical  Motor  Ignition  Traces  (from  D-Series  Testing) 

4.  Motor  Y1  Chamber  Pressure  History 

5.  Still  Photographs  of  Motor  Y1  showing  CCB  consumption.  Here,  the  CCB  is  located  on  the  right  hand 
side  of  the  image  and  flow  is  from  right  to  left.  From  top  to  bottom,  images  are  for  t  =  3.3,  5.45,  5.5, 
and  11.5  sec 

6.  Measured  Chamber  Pressure  Histories  for  Motors  F4,  H3,  and  Ml  at  Average  Pressures  of  Roughly 
100,  200,  and  400  psi,  Respectively 

7.  LDPE  Fuel  Regression  Data  Showing  Effect  of  Pressure 

8.  Comparison  of  Previous  HP  Investigators  Regression  Rate  Measurements  to  Present  Study 

9.  Effect  of  Aft  Combustion  Length  on  C*  Efficiency  with  Inerts 

10.  Effect  of  L*  on  (7*  Efficiency  with  Inerts 

11.  Effect  of  Chamber  Pressure  on  (7*  Efficiency  with  Inerts 

12.  Waterfall  Plot  of  Chamber  Pressure  from  Motor  F4 

13.  Waterfall  Plot  of  Chamber  Pressure  from  Motor  G3 
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Figure  1:  Cross  Sectional  View  of  the  Combustion  Chamber 


Figure  2:  Test  Facility  Schematic 


Figure  5:  Still  Photographs  of  Motor  Y1  showing  CCB  consumption.  Here,  the  CCB  is  located  on  the  right 
hand  side  of  the  image  and  flow  is  from  right  to  left.  From  top  to  bottom,  images  are  for  t  —  3.3,  5.45,  5.5, 
and  11.5  sec 


69 


Motor  Ml 


Motor  H3 


Motor  F4 


Time  (sec) 


Figure  6:  Measured  Chamber  Pressure  Histories  for  Motors  F4,  H3,  and  Ml  at  Average  Pressures  of  Roughly 
100,  200,  and  400  psi,  Respectively 
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Figure  7:  LDPE  Fuel  Regression  Data  Showing  Effect  of  Pressure 


Mass  Flux,  G  (Ibm/fliV^-s)) 


Figure  8 :  Comparison  of  Previous  HP  Investigators  Regression  Rate  Measurements  to  Present  Study 
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Figure  9:  Effect  of  Aft  Combustion  Length  on  C *  Efficiency  with  Inerts 
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8  Appendix  D  -  3-D  Boundary  Element  Model  Development 


Excerpts  from  Chao,  C.,  Ph.D.  Thesis,  “Three-Dimensional  Nonlinear  Jet  Atomization 
Model,”  Purdue  University,  December,  1998. 
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4.  THREE  DIMENSIONAL  BOUNDARY  ELEMENT  METHOD  FOR  NONLINEAR 

DROPLET  SIMULATION 

4.1  Introduction 

A  three  dimensional  Laplace’s  solver  has  been  developed  in  the  previous  chapter.  The 
implementation  of  this  3-D  solver  for  a  nonlinear  droplet  simulation  will  be  described  in 
this  chapter.  Several  auxiliary  modules  had  to  be  built  to  support  these  simulations. 

First,  a  grid  generation  module  had  to  be  developed  for  the  droplet  geometry.  Grid 
generation  will  define  the  physical  domain  for  the  case  we  will  study.  It  also  specifies  the 
proper  initial  conditions  for  the  each  node  on  the  grid.  Note  that  a  given  value  of  velocity 
potential  {<f>)  or  of  the  normal  derivative  of  velocity  potential  (q)  is  required  for  each  node 
in  order  to  form  a  well  posed  boundary  value  problem.  The  grid  generation  process  is 
described  in  Section  4.2. 

Secondly,  a  procedure  to  calculate  the  normal  vector  ,  surface  curvature  and  tangential 
velocity  at  each  node  is  required.  At  every  node,  there  exists  an  unique  normal  vector 
while  an  infinite  number  of  tangential  vectors  exist  at  every  node  location.  The  normal 
vector  serves  as  a  fundamental  factor  for  calculating  other  surface  properties  such  as  surface 
curvature  and  tangential  velocity.  Surface  curvature  is  defined  as  the  rate  of  change  of 
normal  vector,  and  it  was  introduced  into  the  dynamic  boundary  condition  as  defined  in 
Equation  ??  by  the  relationship  with  the  surface  tension  force.  The  experience  gained  from 
2-D  case  shows  that  accurate  methods  have  to  developed  for  calculating  normal  vector  and 
curvature,  especially  for  curvature  which  involves  the  computation  of  the  second  derivative 
of  surface  properties.  Once  the  normal  vectors  are  known,  the  tangential  velocity  can  be 
calculated.  The  combination  of  the  tangential  velocity  and  the  normal  velocity  forms  the 
total  velocity.  The  total  velocity  will  facilitate  the  movement  of  a  node  to  a  new  location. 
The  procedure  to  calculate  normal  vectors,  curvatures  and  tangential  velocities  for  each 
node  are  explained  in  Section  4.3.1. 

The  chapter  is  concluded  with  a  grid  function  convergence  test  (Section  4.4)  and  sample 
simulations  (Section  4.5) 

4.2  Grid  Generation  for  a  Sphere 

The  best  way  to  construct  grids  for  a  sphere  is  to  use  a  polyhedra  with  high  degree 
of  symmetry  and  order  as  a  basic  model.  There  are  five  polyhedra  which  can  be  selected 
as  the  basic  model:  tetrahedron,  octahedron,  cube,  icosahedron,  and  dodecahedron  [1]. 
These  objects  are  all  convex  polyhedrons  and  have  faces  which  are  nonintersecting  regular, 
plane,  convex  polygons  with  straight  edges.  An  equal  number  of  congruent  faces  meet 
at  each  vertex  of  a  particular  polygon,  and  each  has  a  circumsphere  which  touches  all  of 
its  vertices.  Also,  each  polygon  has  the  same  size  as  an  equilateral  triangular  face.  The 
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octahedron  has  been  chosen  as  basic  model  due  to  the  ability  to  easily  find  the  positions  of 
its  six  vertices  which  simply  lie  on  the  x,  y.  and  z  axis  with  the  length  of  radius  in  both 
positive  and  negative  direction  of  each  axis. 

Based  on  this  polygon,  we  are  able  to  construct  finer  grids  by  subdivision  and  projec¬ 
tion.  Subdivision  is  used  to  break  down  an  existing  polygon  in  an  orderly  manner  with  a 
series  of  lines.  The  new  points  produced  by  the  intersection  of  these  lines  will  become  the 
new  vertices.  Figure  4.1  depicts  the  process  of  subdivision  on  one  face  of  an  octahedron. 
The  next  step,  projection,  is  to  imagine  that  the  subdivided  polyhedron  is  surrounded  by 
its  circumscribing  sphere,  with  its  vertices  touching  that  sphere.  Each  edge  of  the  circum¬ 
scribed  polyhedron,  together  with  each  of  the  lines  which  subdivides  its  faces,  should  be 
projected  from  the  center  of  the  polyhedron  onto  the  surface  of  the  circumscribing  sphere, 
as  in  Figure  4.2.  When  the  subdivided  polyhedron  is  projected  onto  the  circumscribing 
sphere,  each  of  its  edges  and  each  line  drawn  across  a  face  becomes  longer.  Some  of  those 
lines  are  projected  further  outwards  and  become  longer  than  the  others,  so  the  edges  of 
the  polyhedron  are  not  all  the  same  length  and  many  triangular  faces  are  not  equilateral. 
Figure  4.3  shows  four  sets  of  spherical  grid  setups. 

4.3  Free  Surface  Module 

There  are  two  methods  to  calculate  the  free  surface  properties.  One  involves  in  prede¬ 
termining  a  given  fitting  function  for  free  surface  and  apply  the  least  square  fitting  scheme 
to  calculate  the  surface  properties.  Another  is  based  on  the  equation  derived  by  [2]  and 
a  suitable  weighting  functions.  The  following  two  subsections  are  devoted  to  these  two 
methods.  The  method  uses  in  3D  model  is  the  first  one  due  to  its  simplicity. 

4.3.1  Calculate  Surface  Properties  by  a  Given  Function 

The  surface  module  will  provide  the  calculation  of  normal  vectors,  surface  curvatures, 
and  tangential  velocities.  We  can  calculate  those  quantities  using  a  fitting  scheme.  The 
calculation  of  the  normal  vector,  the  local  curvature,  and  the  local  tangential  velocity  can 
be  obtained  by  fitting  the  surface  locally  to  a  given  function.  The  least  square  surface 
fitting  will  be  used  as  the  fitting  scheme.  The  following  method  is  following  the  procedure 
by  Chahine  [3]  and  Pozrikidis  [2]. 

Supposed  the  surface  can  be  represented  locally  around  a  given  node  by  a  function 
F(x,y,z)  which  obeys  F(x,y,z )  =  0.  The  second  degree  polynomial  fit  function  is  chosen 
for  the  function  F(x ,  y,  z )  which  has  a  general  form  as: 

F(x,y1z)  =  ax 2  +  by2  +  cz 2  +  fxy  +  gxz  +  hyz  +  px  +  qy  +  rz  +  d  [4.1] 

The  function  F(x ,  y,  z)  is  selected  to  be  the  best  description  the  shape  of  domain  we  are 
interested.  For  example,  for  a  sphere  the  coefficient  d  can  be  set  to  —1  and  /,  g,  and  h  can 
set  to  be  zero.  The  coefficients  p ,  q,  and  r  remain  for  the  case  of  distorted  sphere.  Thus, 
the  fit  function  for  a  sphere  can  be  rewritten  as: 

F(x,  y,  z)  =  aix2  +  a2J/2  +  «3^2  +  +  a$y  +  a6z  -  1  [4.2] 

The  coefficients  cq  to  a6  can  be  obtained  by  direct  fit  method  or  least  square  method. 
For  the  direct  fit  method,  the  coefficients  eq  to  a6  are  determined  by  choosing  six  points 


77 


which  surround  a  given  node.  Numerical  experiments  show  that  this  method  produced  less 
accurate  results.  Thus,  the  direct  fit  method  was  not  used  in  the  simulations.  For  the 
least  square  method,  the  errors  to  fit  the  function  can  be  reduced  to  some  level  by  using 
more  surrounding  nodes  to  calculate  the  coefficients  aj  to  a^.  Depending  on  the  grid,  the 
surrounding  nodes  can  be  used  up  to  the  third  loop  around  a  given  node.  For  example,  for 
the  grid  of  258  nodes,  total  surrounding  nodes  of  a  given  node  which  includes  the  three  loop 
nodes  around  it  will  has  about  30  nodes.  The  coefficients  a\  to  a6  are  then  determined  by 
a  least  square  fit  to  the  surface  for  these  30  nodes.  Once  the  fitting  function,  F{x ,  y,  z),  is 
obtained,  the  local  normal  vector  can  be  calculated  by  : 


n  =  ± 


VF 

|VF| 


[4.3] 


The  appropriate  sign  is  chosen  to  insure  the  normal  is  pointed  outward.  Writing  the  normal 
vector  by: 

n  =  mi  +  n2j  +  n3fc  [4.4] 


then: 


Ul  ~  |VF| ’  n2-|VF|’  ns_|VF| 


[4.5] 


where  Fx ,  Fy  and  Fz  are  the  first  derivative  of  function  F(x,y,z)  with  respect  to  a:,  y, 
and  z,  respectively.  |VF|  is  the  norm  of  gradient  of  function  F(x,  y,  z)  which  is  defined  as 

t  . -  /V  A  A 

Jn\  +  n\  +  n§.  And  ni,  n2  and  n3  are  three  components  of  normal  vector  n  in  i,  j,  and  k 
direction,  respectively.  The  local  surface  curvature  is  obtained  by: 


k  =  V  •  n 


[4-6] 


Introducing  Equation  4.2,  curvature  is  given  by: 

FrAFl  +  FT)  Fw(Fj  +  F?)  F„(i*  +  F;> 

|VF|3  T  |VF|3  |VF|3 

where  Fxx,  Fyy,  and  Fzz  are  the  second  derivative  of  function  F(x,y,z)  with  respect  to  x, 
y,  and  z,  respectively.  In  order  to  calculate  the  tangential  velocity,  the  fitting  process  is 
performed  on  the  velocity  potential  in  a  similar  manner  as  for  the  normal  vector.  Since 
the  velocity  potentials  are  given  on  each  vertex,  the  process  of  the  surface  fitting  on  <j> 
is  analagous  to  our  approach  for  k.  The  polynomial  function  has  the  similar  form  as  for 
function  F(x ,  y,  z)  and  we  denotes  as  G(x,  y,  z)  and  is  given  by: 

( f>{x ,  y,  z)  «  G(x ,  y,  z)  =  bix 2  +  62y2  +  b3z2  +  b4x  +  b5y  +  b6z  [4.8] 

The  coefficients  b\  to  b$  are  determined  by  least  square  method.  Once  the  function  <j>(x,  y,  z) 
is  determined  locally,  the  tangential  velocity  can  be  calculated  by: 

vt  =  h  x  (V<£  X  h)  [4.9] 

After  introduction  Equation  4.2  and  Equation  4.8,  tangential  velocity  can  be  obtained  by: 

vt  =  vn  i  +  vt2j  +  vak  [4.10] 


78 


with: 


where 


Vtl  = 

^GnS  ~  ^3^712 

[4.11] 

Vt2  = 

^3Gn2  —  n\Gn3 

[4.12] 

Vt3  = 

niGn2  —  U2Gn\ 

[4.13] 

Gn  i  = 

Gyn$  —  Gzri2 

[4.14] 

Gn2  = 

G  zTl\  G  3 

[4.15] 

GnZ  — 

Gx  ft  2  ~  Gy  Til 

[4.16] 

Gx,  Gy,  and  Gn  are  the  first  derivative  of  the  function  G(x,y,z)  with  respect  to  x,  y,  and 
x,  respectively.  Thus,  by  combing  the  tangential  velocity  and  the  normal  velocity,  the  total 
velocity  can  be  determined  and  is  given  by: 

v  =  vt  +  qn  [4-17] 

where  q  is  the  normal  derivative  of  velocity  potential. 


4.3.2  Calculating  Surface  Properties  Using  Weighting  Functions 

A  more  general  method  to  calculate  surface  properties  not  based  on  a  predetermined 
function  is  described  by  Pozrikidis  [2].  The  basic  idea  involves  using  the  weighting  average 
methods  to  calculate  surface  properties.  Normal  vectors  for  each  element  are  calculated  first. 
Since  normal  vectors  of  each  triangular  element  can  be  calculated  easily  by  the  information 
given  by  the  three  vertices’  position.  The  normal  vectors  for  each  node  are  taking  a  weighted 
ng  average  over  the  surrounding  elements  about  each  node  which  is  given  by: 


e 


*4  =  £ 

2=1 

e 

[4.18] 

ny  —  ny^ 

t= 1 

e 

[4.19] 

»i  =  £<n’ 

i= 1 

[4.20] 

where  nlx ,  nx,  and  n\ j,  are  the  three  components  of  normal  vector  a  node  i  and  ST  is  the 
weighting  function.  Here,  e  denotes  the  total  number  of  element  surrounding  node  i.  The 
success  of  this  method  relies  on  the  selection  of  a  suitable  weighting  function.  The  first  try 
for  selection  of  the  weighting  function  is  to  use  the  arithmetic  averaging.  The  implementa¬ 
tion  of  arithmetic  averaging  will  result  in  twice  averaging  the  normal  vectors  which  could 
introduce  too  much  damping.  The  numerical  results  proved  this  speculation  by  reducing 
droplet  movement  when  compared  to  the  technique  described  in  Section  4.3.1.  Thus,  a 
more  advanced  weighting  method  was  desired. 

The  curvature  can  be  calculated  by  the  method  described  by  Pozrikidis  [2].  The  mean 
curvature  (km)  at  a  given  point  is  given  by: 

km  =  |(*('  0)  +  k{oo)) 


[4.21] 
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*(0)  = 


1  dn  dx 
a™  dr]  dr] 


.  .  .  1  dn  dx 

k{  oo)  =  — -7TT- 
avv  dr]  dr] 


[4.23] 


where  the  metric  tensor  a  is  given  by: 


a«  = 


dx  dx 
dr]  dr] 
dx  dx 
dr]  d£ 
dx  dx 

d£  di 


[4.25] 

[4.26] 


And,  by  using  linear  triangular  element,  the  surface  properties  are  given  by: 

P(v,  0  =  (1  -  V  -  QPi  +  vP2  +  £ft  [4.27] 

Thus,  the  metric  tensor  is  given  by: 

o-m  =  xl  +  Vv  +  zv  t4-28] 

a-nt  =  +  ylv\  +  z\z\  [4.29] 

a(,i  =  x\  "h  2/|  +  2|  [4.30] 

Substituting  the  above  three  equations  into  Equation  4.22  and  Equation  4.23,  we  obtain; 


*(0)  = 

fc(oo)  = 


CL-rin  (^£,77  *<)  +  4"  nz,r\zri) 

1 


[4.31] 

[4.32] 


The  mean  curvature  is  calculated  by  Equation  4.21. 

The  velocity  potential  is  given  at  each  node  and  its  distribution  over  an  element  is  given 
by: 

<£(»?, 6  =  (1  “  V-0<t>i  +  V<h  +  €<h  [4-33] 

Thus,  the  components  of  the  gradient  potential  velocity  for  a  given  element  are  given  by: 


d<f>  _  d<p  dr]  d(j>  d£ 
dx  dr]  dx  ^  d£  dx 
d4>  _  d(f>  dr]  d(f>  d£ 
dy  dr]  dy  +  di  dy 
d<p  _  d<j)  dr]  d(j>  d£ 
dz  dr]  dz  d£  dz 


[4.34] 

[4.35] 

[4.36] 


4>r,  =  <h  ~  <fri  ;  fa  =  -  <t>\  [4-37] 

In  order  to  obtain  those  components  of  the  gradient  potential  velocity  at  each  node,  an 
averaging  scheme  needed  to  be  implemented.  Once  the  gradient  potential  at  each  node  is 
calculated,  the  the  tangential  velocity  and  total  velocity  can  be  calculated  by  Equation  4.17 
and  Equation  4.9,  respectively. 
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Table  4.1  Grid  Setup  for  the  Test  of  Grid  Convergence 


Grid  Set 

No.  of  Nodes 

No.  of  Elements 

No.  1 

6 

8 

No.  2 

18 

32 

No.  3 

66 

128 

No.  4 

258 

512 

4.4  Grid  Convergence  Test 

Since  the  3-D  model  is  developed  to  solve  the  Laplace  equation,  there  exists  an  analytical 
solution  for  a  simple  domain  like  a  sphere  by  specifying  the  boundary  conditions  on  the 
surface  of  the  sphere  of  either  Dirichlet  or  Neumann  type.  For  the  testing  given  below,  the 
boundary  condition  of  the  Dirichlet  type  has  been  used. 

The  boundary  condition  is  given  by  specifying  a  function,  say  $(x,y,  z),  which  will 
satisfy  the  Laplace  equation.  Then,  the  normal  derivatives  of  $(£,  y,  z)  can  be  calculated 
analytically  by  qana  —  d$/dn  and  the  results  will  be  used  to  compare  the  results  from  3-D 
BEM  model.  The  subscript  of  qana  denotes  the  analytical  solution  for  q.  The  function  for 
y)  z)  has  been  chosen  to  be: 

—  exp(x)cos(y)  +  exp(z)sin(x)  [4.38] 

And,  for  a  sphere,  the  surface  can  be  represented  by: 

F{x,  y,z)  =  x2  +  y2  +  z2  -  1  [4.39] 

which  provides  the  qana  can  be  calculated  by: 

9ono(®>  Vi  ■z)  =  =  ^  [4.40] 

The  h  values  are  calculated  by  V77/ 1 V /']  described  in  Section  4.3.1.  Four  sets  of  grids  have 
been  used  for  the  test  of  grid  convergence.  Table  4.1  is  highlights  the  number  of  nodes  and 
elements  used  for  the  grid  convergence  test. 

Figures  4.4,  4.5,  4.6,  and  4.7  show  the  error  in  q  between  numerical  and  analytical 
solutions  for  these  four  different  grids.  For  the  grid  with  only  six  nodes  on  the  spherical 
surface,  the  difference  between  the  analytical  and  numerical  results  shows  a  large  discrep¬ 
ancy  in  Figure  4.4.  However,  the  error  will  reduce  about  one  order  when  a  finer  grid  setup 
with  18  nodes  on  the  spherical  surface  is  used.  The  trend  of  reducing  error  is  held  for  using 
the  even  more  finer  grid  setup  as  66  nodes  and  258  nodes  are  used.  Figure  4.8  is  a  error 
plot  for  these  four  different  grid  setups. 

Another  way  to  look  at  the  grid  convergence  phenomenon  is  by  showing  the  error  re¬ 
ducing  trend  for  choosing  the  same  grid  points  which  appear  on  the  these  four  grid  setups. 
Table  4.2  indicates  that  there  are  six  points  which  appear  in  all  the  four  grids. 

These  six  points  will  be  remained  when  the  subdivision  and  projection  process  be  applied 
to  construct  a  finer  grid  for  a  sphere.  Table  4.3  shows  the  error  decreases  as  the  finer  grids 


Table  4.2  Six  Common  Grid  Points  for  the  our  Different  Grid  Setups 


Node  Location 

X 

y 

z 

No.l 

0.0 

1.0 

0.0 

No.2 

0.0 

0.0 

1.0 

No. 3 

1.0 

0.0 

0.0 

No.4 

-1.0 

0.0 

0.0 

No.5 

0.0 

-1.0 

0.0 

No.6 

0.0 

0.0 

-1.0 

Table  4.3  Results  of  Test  of  Grid  Convergence:Part  I 


qerr 

Node  Location 
No.l 

Node  Location 
No.2 

Node  Location 
No.3 

6  nodes 

0.266241 

-0.032797 

-0.864719 

18  nodes 

-0.032538 

-0.021071 

-0.037991 

66  nodes 

-0.005458 

-0.002019 

0.055933 

258  nodes 

-0.003276 

-0.000160 

0.024899 

qerr 

Node  Location 
No.4 

Node  Location 
No.5 

Node  Location 
No.6 

6  nodes 

-0.269631 

0.266241 

-0.032797 

18  nodes 

0.050957 

-0.032528 

-0.021071 

66  nodes 

0.018494 

-0.005458 

-0.002019 

258  nodes 

0.005620 

-0.003276 

-0.000160 
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are  applied  for  all  of  the  six  common  points.  While  the  scheme  is  formally  second  order, 
the  convergence  rates  are  not  as  rapid  as  one  would  expect.  Here,  the  largest  error  result 
be  attributed  to  the  changes  in  surface  shape  with  changing  number  of  nodes.  In  this  case, 
convergence  is  also  measuring  the  rate  at  which  the  grid  converges  to  a  perfect  sphere. 
Based  on  the  evidence  in  these  calculations  (and  in  Chapter  ??)  the  code  was  presumed  to 
be  valid  for  use  in  time-dependent  simulations. 


4.5  Sample  Calculation  -  Droplet 


In  order  to  verify  the  time-dependent  capabilities  of  the  3-D  code  developed  in  the  pre¬ 
vious  sections,  the  model  was  used  to  simulate  the  distortion  of  a  droplet.  The  driving 
force  for  the  distortion  of  a  droplet  is  due  to  surface  tension.  Different  modes  oscillation 
are  investigated  for  the  droplet  simulation.  In  order  to  show  the  different  modes  movement 
of  droplet,  the  initial  values  of  velocity  potential  are  specified  by  Legendre  function.  The 
Legendre  function  is  the  solution  of  Laplace  equation  in  spherical  polar  coordinate  system 
and  symmetry  about  the  z  axis,  ie.  the  dependent  variable  of  Laplace  equation  is  indepen¬ 
dent  of  6  where  6  £  [0,2tt].  Thus,  the  Legendre  functions  are  the  solution  of  the  following 


equation: 


V2u 


d2u  2du  1  d2u  cot  <j>f  du 
dr 2  r  dr  r2  d<f>f 2  r2  d(j>f 


[4.41] 


where  <j/  £  [0,  27t)  and  the  prime  added  to  distinguish  the  variable  used  for  the  velocity 
potential  <f>.  The  dependent  variable  is  velocity  potential  in  this  case.  Therefore,  the  initial 
value  for  velocity  potentials  are  specified  by: 


<£(r,  0,  <j>r)  =  ePn(cosO) 


[4.42] 


with, 

0<(j><7r  0  <  0  <  2tt  [4.43] 

where  Pn(cos0)  is  the  Legendre  function  of  order  n .  The  order  will  refer  to  the  mode 
number  for  the  following  simulation. 

Figure  4.9  shows  the  half  period  of  third  mode  oscillation  for  e  =  0.001  with  66  nodes 
and  128  elements  on  the  spherical  surface.  The  characteristic  of  the  third  mode  oscillation 
is  the  three  lobe  shape  on  the  spherical  surface.  Figure  4.10  shows  the  half  period  of  the 
forth  mode  oscillation  with  the  same  conditions  as  the  third  mode  oscillation.  The  four  lobe 
shape  is  developed  for  the  fourth  mode  oscillation.  Since  no  regridding  process  has  been 
applied  in  these  cases,  the  simulations  were  restricted  to  a  small  perturbations. 


Subdivided  by  two 


Figure  4.1  Subdivision  of  a  Face  of  an  Octahedron 
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Figure  4.2  Projection  of  Subdivided  Nodes  on  the  Surface  of  a  Sphere 
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.5  Compariso 


>n  of  Analytic  and  Numerical  Surface  Velocities  for  a  18  Node  G 
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total  node  =  258,  total  elements  =  512 


node  number 


Figure  4.7  Comparison  of  Analytic  and  Numerical  Surface  Velocities  for  a  258  Node  Grid 
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Figure  4.8  Standard  Deviation  in  q  vs.  Four  Different  Grid  Setups:  6,  18,  66,  and  258 
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5.  THREE  DIMENSIONAL  BOUNDARY  ELEMENT  METHOD  FOR  NONLINEAR 

INFINITE  LIQUID  JET  SIMULATION 

5.1  Introduction 

A  process  to  solve  the  problem  of  nonlinear  deformation  of  an  infinite  jet  subjected 
to  both  axial  and  circumferential  disturbances  by  using  3-D  BEM  model  developed  in 
Chapter  ??  is  described  in  this  chapter.  This  process  begins  with  a  description  of  an 
analytic  temporal  instability  analysis  by  Yang  [4]  (as  described  in  Section  5.2).  The  resulting 
dispersion  equation  has  been  simplified  to  fit  our  assumptions  and  serves  as  an  analytic 
solution  to  compare  to  the  numerical  results  from  3-D  BEM  model.  Section  5.3  describes 
a  grid  generation  process  for  the  liquid  column.  The  liquid  column  serves  as  a  portion  of 
an  infinite  liquid  jet.  The  initial  shape  of  jet  can  be  arbitrary  specified  by  manipulated  the 
longitudinal  (streamwise)  wave  number  and  transverse  (circumferential)  wave  numbers. 

As  in  the  previous  chapter,  a  free  surface  module  is  developed  specifically  for  the  jet 
geometry  in  Section  5.4.  Different  fit  functions  are  chosen  for  the  jet  due  to  the  differ¬ 
ent  geometric  shape.  The  process  of  dealing  with  the  corner  problem  is  described  in  this 
section.  A  view  of  nonaxisymmetric  wave  on  the  liquid  jet  is  described  in  Section  5.5.  Sec¬ 
tion  5.6  presents  the  results  of  the  comparison  is  made  between  the  analytic  linear  stability 
analysis  (Section  5.2)  numerical  model. 

5.2  Dispersion  Equation 

The  mathematical  model  used  for  the  temporal  instability  analysis  is  based  on  linear 
analysis  conducted  by  Yang  [4]. The  resulting  dispersion  equation  has  been  simplified  to 
fit  our  assumptions  of  inviscid  flow  with  negligible  gas-phase  effects.  Surface  tension  is  a 
dominant  force  for  the  instability  of  the  low  speed  jet  which  is  also  called  the  capillary  jet. 
Rayleigh  [5]  showed  that  only  unstable  mechanism  to  cause  the  breakup  of  a  low  speed  jet 
is  the  growth  of  axisymmetric  waves  which  are  also  referred  to  as  dilation  waves.  For  the 
high  speed  jet,  the  transverse  modes  are  introduced  by  viscous  and  turbulent  effects  at  the 
nozzle  exit.  Asymmetric  waves  are  also  called  sinuous  waves.  In  order  to  understand  the 
behavior  of  the  sinuous  wave,  the  theory  has  to  include  the  transverse  mode  oscillation. 

Yang’s  theory  begins  with  a  two  stream  flow,  liquid  jet  and  gas  stream,  which  are 
assumed  to  be  incompressible  and  inviscid  fluid.  Liquid  jet  with  density  p\  and  issued  at  a 
uniform  speed  U\,  while  a  coaxial  gas  with  density  P2  flowed  with  a  uniform  speed  U2.  the 
interface  between  liquid  and  gas  phase  has  been  perturbed  by: 

v  =  m  =  m  =  noei{kz+mB)+at  [5.1] 

where  770  is  perturbation  magnitude,  k  and  m  are  the  wave  number  in  streamwise  or  longi¬ 
tudinal  and  azimuthal  or  transversal  direction,  respectively.  Here,  a  is  a  complex  variable, 
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a  =  o.r  +  ia%.  The  real  part  of  a,  ar,  is  the  growth  rate.  The  positive  valve  of  ar  indicates 
the  wave  growth,  while  the  negative  value  indicates  the  wave  is  stable  and  will  eventually 
damp  out.  The  imaginary  part  of  a,  <*;,  is  connected  to  the  wave  propagation  speed.  The 
resulting  characteristic  equation,  i.e.,  dispersion  equation  is  given  by: 

/  *\2  _  .  fca  1  —  TO  —  (ka)  r,,  21 

[0lr)m  ~  (7 m  +  PmQ )2  We  7m  +  PmQ 


where: 


K)m  = 
We  = 
Q  = 

T m  = 

Pm  — 


(Qr)^ 

(C/i  -  W/«2 

q(£/i  -  U2fpx) 

a 

P2 

Pi 

klm  (ka) 
I'mika) 
kKwi  (ka) 
Kin(ka) 


[5.3] 

[5.4] 

[5.5] 

[5.6] 

[5.7] 


And,  a  is  the  radius  of  nozzle.  Im  is  the  modified  Bessel  function  of  the  first  kind  and 
order  m.  Km  is  the  modified  Bessel  function  of  the  second  kind  and  order  m.  The  prime 
on  Im  and  Km  denotes  the  first  derivative  with  respective  to  their  parameters.  For  an 
axisymmetric  wave  oscillation,  the  azimuthal  wave  number  m  is  simply  equal  to  zero.  Thus, 
the  asymmetric  instability  is  introduced  by  providing  nonzero  value  of  m. 

In  order  to  fit  in  our  previous  assumptions,  we  let  Q  and  U2  be  equal  to  zero  for  the 
absence  of  gas  phase  which  also  provides  We  =  1.  Thus,  we  have 


(«r )m  =  (Ht1  -  ™2 


(ka2)] 


tN 

klm (ka) 


[5.8] 


Figure  5.3  shows  the  plot  of  growth  rate,  (a*)2m  vs  ka  at  different  values  of  m.  It  indicates 
that  the  only  unstable  wave  for  a  quiescent,  incompressible,  and  inviscid  liquid  jet  is  ax¬ 
isymmetric  wave.  All  of  non-axisymmetric  modes  (m  >  0)  are  stable  due  to  the  negative 
values  of  growth  rate  over  all  of  dimensionless  streamwise  wave  number  ka. 

Thus,  the  dispersion  equation  is  used  for  the  temporal  instability  analysis  later  in  this 
chapter  is  given  by: 

K)S  =  (*«)[! P-91 

which  is  identical  to  the  result  of  Rayleigh’s  analysis. 


5.3  Grid  Generation  of  a  Liquid  Column 

The  first  step  in  the  numerical  analysis  is  to  set  up  the  grid  points  for  a  liquid  jet.  A 
cylinder  with  a  small  perturbation  on  the  surface  is  used  to  represent  the  initial  shape  of 
the  liquid  jet.  Figure  5.1  depicts  the  small  perturbation  on  the  surface  of  the  liquid  jet  in 
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Table  5.1  Three  Different  Grid  Setups  for  a  Liquid  Column 


int_z 

int.the 

int_r 

total  nodes 

total  elements 

Set  1  5 

5 

2 

42 

80 

Set  2  20 

6 

3 

222 

440 

Set  3  30 

10 

3 

452 

900 

the  r  -  z  plane  and  r  -  9  plane.  Since  the  liquid  jet  is  assumed  to  be  of  infinite  length,  the 
computational  domain  is  chosen  to  be  one  wavelength  (A)  of  the  perturbation. 

A  grid  generation  scheme  has  been  applied  on  a  portion  of  an  infinite  liquid  jet  which  is 
described  by  a  cylindrical  column.  This  cylindrical  column  is  constituted  by  a  side  surface 
and  two  end  faces.  Each  surface  requires  a  different  method  to  construct  grid  points  and 
elements.  The  result  of  grid  generation  is  to  provide  information  about  the  location  of  each 
grid  points  on  the  surface  and  what  three  nodes  to  make  up  each  triangular  element. 

The  number  of  intervals  in  the  r,  0 ,  and  z  directions  (int_r,  int_z,  and  int_the,  respec¬ 
tively)  must  be  input  in  order  to  generate  the  grid.  The  value  of  int_z  will  show  int_z  + 
1  nodes  on  the  wave  along  the  longitudinal  direction,  while  int.the  +  1  nodes  along  the 
transversel  direction.  The  total  length  of  a  liquid  column  is  given  by  n  •  A  where  n  is  the 
number  of  wavelengths  and  A  denotes  wavelength,  note  that  A  is  given  by  A  =  2ir/k  where 
k  is  the  wave  number  in  the  streamwise  or  longitudinal  direction.  The  initial  shape  of  the 
liquid  column  is  given  by: 

r  =  1  +  8cos{kz  +  md )  [5.10] 

where  8  is  the  perturbation  magnitude.  For  an  undisturbed  column,  5  is  given  by  zero, 
and  an  axisymmetric  column  is  given  by  m  =  0.  Table  5.3  highlights  three  grid  setups,  and 
the  geometrical  perspective  of  there  three  grid  setups  are  shown  on  Figure  5.4. 

5.4  Free  Surface  Modules 

The  process  to  calculate  the  surface  normal  vector,  surface  curvature,  and  tangential 
velocity  is  the  same  as  the  methods  used  in  the  droplet  case  as  discussed  in  Chapter  4. 
There  are  three  major  differences  between  the  droplet  case  and  the  jet  case.  First  is  the 
usage  of  a  different  fit  function  for  the  surface  fitting  due  to  the  different  geometrical  domain 
in  the  liquid  jet  case.  The  least  square  method  is  also  applied  for  the  fitting  scheme  in  the 
liquid  jet  case,  but  here  the  fit  function  is  chosen  to  be: 

F(x,  y,  z)  =  a\x2  +  a2y2  +  a3x  +  a4y  +  a5z  -  1  [5.11] 

And,  the  fit  function  for  the  velocity  potential  fitting  in  the  liquid  jet  case  is  chosen  to  be: 

<t>{x,  y ,  z)  =  G(x,  y,  z)  =  b\x2  +  b2y 2  +  b3x  +  b4y  +  b5z  [5.12] 

The  accuracy  of  the  fitting  function  can  be  determined  by  observing  the  results  of  surface 
curvature  for  an  undisturbed  jet.  For  an  undisturbed  jet  which  is  a  smooth  cylindrical 
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column,  the  surface  curvature  is  equal  to  one  at  every  node  on  the  surface  side  which  does 
not  include  the  nodes  on  the  two  end  faces. 

The  second  difference  between  droplet  case  and  jet  case  is  the  appearance  of  the  corner 
nodes  in  the  jet  case.  The  corner  nodes  are  referred  to  the  nodes  that  lie  on  the  intersection 
potion  between  the  surface  side  and  two  end  faces.  Those  nodes  have  multiple  value  of 
q  which  is  q  —  d(f>/dn.  Since  these  two  end  faces  are  the  cross  sections  of  an  infinite  jet 
which  serve  as  symmetry  planes.  By  the  definition,  the  values  of  q  are  set  to  be  zero  for  the 
nodes  on  these  symmetry  planes.  Thus,  the  elements  lying  on  one  of  the  two  end  faces  will 
have  no  contribution  to  the  multiple  value  problem  of  q  since  q  is  zero  on  these  elements. 
However,  for  the  elements  on  the  surface  side,  an  assumption  is  made  to  deal  with  the 
multiple  value  of  q.  Since  the  multiple  value  problem  is  arisen  due  to  the  different  normal 
vectors  associated  with  those  surrounding  elements,  these  differences  in  the  normal  vectors 
can  be  assumed  to  be  small  when  a  finer  grid  is  used.  Therefore,  multiple  value  of  q  for 
those  elements  becomes  a  single  value  of  q. 

An  example  is  given  to  further  explain  how  to  deal  with  the  multiple  value  of  q  on 
the  corner  nodes.  For  grid  setup  number  1  mentioned  in  Table  5.3,  the  corner  node  and 
their  surrounding  elements  are  shown  on  the  Figure  5.5.  The  zoom-in  picture  in  Figure  5.5 
shows  the  node  numbers  and  element  numbers  surrounding  node  6.  The  dark  portion  of 
this  zoom-in  picture  indicates  the  end  face  while  the  white  portion  of  this  picture  refers 
to  the  surface  side  of  column.  The  summation  process  in  Equation  ??  gives  the  following 
equation  around  node  6: 

...  =  - 1-  [910.3*^3, 6.10  +  949.3*^3,6.49  +  950.2*^2,6.50 

+967.1*^1,6  .67  -  967.1*^2,6.67  ~  967.1 *^3, 6. 67 
+968.1*^1,6  .68  ~  968.1*^2, 6.68  ~  968.1*^3,6.68 

+980.2*^2,6.80  +  ’  '  •]  [5.13] 

The  notation  of  qlj  k  is  that  i  refers  to  node  number,  j  indicates  the  element  number  and  k 
is  the  node  order  in  the  j  element.  For  example,  gfs.i  indicates  the  q  value  at  node  6  which 
lies  on  the  first  node  of  element  68.  And,  is  the  value  of  kernel  function  Si  which  the 

integration  is  performed  on  the  base  point  m  over  element  n.  Thus,  £3,6.10  is  the  value  of 
kernel  function  S3  which  is  based  on  the  base  point  6  and  integrating  over  element  10.  In 
this  example  as  shown  on  Figure  5.5,  we  have: 

•  node  12,  6,  and  30  are  corner  nodes. 

•  node  5  and  29  are  the  nodes  on  the  surface  side  of  the  jet. 

•  node  39  and  38  are  the  nodes  on  the  end  face  of  the  jet. 

•  element  10,  49,  and  50  lies  on  the  surface  side  of  the  jet. 

•  element  68,  67  and  80  lie  on  the  end  face  of  the  jet 
According  to  our  discussion  and  assumption,  we  have:  By  assumption: 


9io.3  —  950.2  —  949.3  —  96 


[5.14] 
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By  definition: 

<767.1  =  <768.1  =  <7so.i  =  0  [5.15] 

Thus,  Equation  5.13  becomes: 

•  .  .  - - 1-  g6[S3,6.10  +  53,6.49  +  52,6.5o]  H -  [5.16] 

The  last  difference  between  the  droplet  case  and  jet  case  is  that  special  attention  is 
warranted  on  the  corner  nodes  for  the  fitting  process.  The  adjacent  nodes  around  a  corner 
node,  which  are  used  for  the  least  square  method,  can  not  be  used  if  they  lie  on  the  end  face 
of  the  jet  due  to  the  dramatic  change  on  the  geometrical  properties  such  as  normal  vector 
and  curvature.  In  this  case,  the  image  points  associated  with  the  nodes  on  the  surface  side 
of  the  jet  which  are  around  a  corner  node  are  used.  The  result  shown  for  the  curvature  for 
an  undisturbed  jet,  the  value  of  one  are  given  on  every  corner  nodes  by  using  the  image 
point  method. 

5.5  Nonaxisymmetric  Modes 

The  initial  shape  of  a  liquid  column  is  given  by  Equation  5.10.  By  varying  the  value  of 
longitudinal  wave  number(A;),  transverse  wave  number(m),  and  perturbation  magnitude^), 
one  can  see  some  asymmetric  shapes  due  to  the  presence  of  the  nonaxisymmetric  mode. 

Figure  5.6  shows  the  column  shape  for  different  m  values  at  a  given  k  —  1.  The  total 
column  length  is  three  wavelengths  and  S  =  0.2.  An  axisymmetric  shape  is  shown  on  the 
first  picture  of  Figure  5.6  which  m  —  0.  The  radius  of  circular  shape  of  cross  section  varies 
along  with  the  longitudinal  direction.  This  type  of  axisymmetric  variation  shape  is  called 
dilation  or  varicose  mode.  For  the  next  picture  of  Figure  5.6,  m  =  1,  the  cross  section  of 
column  is  nearly  circle  and  the  size  of  cross  section  is  constant  along  with  the  longitudinal 
direction.  This  shape  is  called  sinuous  or  snake  mode.  It  shows  that  the  asymmetric  mode 
is  taking  its  effect  on  the  formation  of  column  shape.  The  cross  section  becomes  elliptic  for 
m  =  2  mode  on  the  third  picture  of  Figure  5.6.  And  the  m  lobes  shape  will  appears  on  the 
cross  section  when  different  values  of  m  been  applied. 

For  the  m  =  1  mode  on  the  Figure  5.6,  there  is  an  interesting  shape  variation  when 
the  perturbation  magnitude  varies.  Figure  5.7  shows  four  cross  section  pictures  for  four 
different  perturbation  magnitudes.  The  original  nearly  circular  shape  at  S  =  0.2  becomes 
to  squeeze  at  the  bottom  of  nearly  circular  shape  when  6  increases.  This  clearly  shows  an 
example  of  non-axisymmetric  deformation. 

Figure  5.8  shows  the  cross  section  and  column  shape  for  m  =  2  mode  when  perturbation 
magnitude  is  increased.  The  development  of  two  lobe  shapes  from  the  cross  section  becomes 
clearer  when  S  is  increased.  At  large  5,  the  column  shape  just  like  a  two  cylindrical  columns 
twisted  together. 

Figure  5.9  show  the  cross  section  and  column  shape  for  m  =  3  mode  at  increasing  per¬ 
turbation  magnitude.  The  three  lobe  shape  shown  on  the  cross  Section  as  the  characteristic 
of  the  third  transverse  mode  oscillation.  At  large  5 ,  the  column  shape  looks  like  three 
cylindrical  columns  twisted  together. 
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5.6  Results  of  Linear  Analysis 

The  mathematical  model  developed  by  Yang  and  described  in  section  5.2,  is  used  for  the 
temporal  instability  analysis  on  a  liquid  jet  which  is  assumed  as  incompressible,  inviscid, 
and  quiescent  fluid.  The  dispersion  Equation  5.9  serves  as  the  analytical  solution  for  the 
growth  rate.  The  analytical  solution  is  used  to  compare  to  the  results  from  3D  BEM 
model.  The  calculation  of  growth  rate  for  the  numerical  result  is  based  on  the  derivation 
by  Mansour  [6]  which  he  utilized  the  prediction  of  Rayleigh’s  analysis.  The  extension  of 
derivation  by  Mansour  for  3D  jet  is  described  as  follow. 

The  development  of  jet  is  assumed  to  be  in  the  form  of: 

r  =  1  +  8cos(kz  +  m6)cosh(art )  [5-17] 


where  ay  is  the  growth  rate  which  is  the  real  part  of  a  which  is  defined  previously.  And, 
the  velocity  in  radius  direction  (14)  is  calculated  by: 


—  =  8cos(kz  +  m6)sinh(art)ar 


[5.18] 


Thus  the  growth  rate  ay  is  given  by: 


(r  —  l)2  —  82cos2(kz  +  mO) 


[5.19] 


Since  the  only  unstable  wave  in  the  incompressible,  inviscid,  and  quiescent  fluid  is  the 
axisymmetric  wave,  m  is  set  to  be  zero.  Along  the  longitudinal  direction,  the  nodes  lying 
on  the  trough  of  wave  are  chosen  to  be  the  location  which  the  temporal  analysis  is  conducted 
(thus  kz  =  7r).  Therefore,  the  preceding  equation  is  simplified  to  the  equation  which  was 
used  in  the  Mansour’s  analysis  which  is  given  by: 


o  V2 

a2  = - — - 

r  (r  -  l)2  -  S2 


[5.20] 


In  this  case, only  the  axisymmetric  wave  develops,  the  nodes  on  the  surface  side  of  jet  will 
move  only  in  the  radial  direction  and  Vr  is  given  by  q  which  is  calculated  by  3D  BEM  model. 
Thus,  the  temporal  instability  analysis  for  axisymmetric  oscillation  can  be  conducted  based 
on  Equation  5.2  for  analytical  solution  and  Equation  5.20  for  numerical  solution. 

The  temporal  instability  analysis  is  performed  by  picking  several  longitudinal  wave  num¬ 
bers  ( k )  and  letting  transverse  wave  number  ( m )  to  b  zero  for  axisymmetric  case.  The  total 
wavelength  is  set  to  be  one.  A  single  grid  setup  is  not  appropriate  to  run  this  analysis  due 
to  the  aspect  ratio  of  grid  element  varying  for  different  values  of  k.  Figure  5.10  shows  the 
grid  setup  with  intjz  =  40,  int.the  =  10,  and  intr  =  3  for  three  different  k  values:  k  =  0.3, 
k  =  0.7  and  k  =  0.9.  The  aspect  ration  AO/Az  will  increases  when  Az  decreases  due  to 
increasing  k  values.  Note  that  A  =;  2k /k  and  Az  =  X/intjz.  The  total  length  of  jet  is  chose 
to  be  one  which  results  in  a  large  aspect  ratio  for  higher  k  value.  The  shape  of  jet  will 
looks  bigger  for  those  higher  k  values  due  to  the  higher  aspect  ratio.  A  grid  setup  with 
nonuniform  aspect  ratio  for  element  will  effects  the  accuracy  of  results.  In  order  to  get 
better  results,  the  optimized  grid  needs  to  be  used  for  each  k  value. 

Figure  5.11  shows  that  the  growth  rate  calculated  by  3D  BEM  model  reaches  the  growth 
rate  predicted  by  Equation  5.9  when  grid  setup  approaches  the  optimized  grid  for  k  =  0.7. 
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The  growth  rate  converges  quickly  to  a  constant.  The  same  behavior  will  show  for  different 
k  values.  The  convergent  constant  will  be  the  desired  value  of  growth  rate. 

Figure  5.12  shows  the  same  plot  for  k  =  0.9.  We  see  that  the  track  of  growth  rate 
vs.  time  comes  down  to  the  predicted  value  when  a  finer  grid  is  used.  The  finest  grid 
been  tested  (int.z  =  70,  intJthe  =  10,  int.r  -  3)  results  in  752  nodes  and  1500  elements 
on  the  surface.  The  even  finer  grid  can  not  be  used  in  our  computer  facility  due  to  the 
shortage  of  memory.  However,  the  convergence  behavior  is  shown  on  Figure  5.12  when  the 
fine  grid  is  implemented.  We  may  expect  that  the  computed  results  will  finally  agree  with 
the  predicted  values  when  the  very  fine  grid  is  used.  Figure  5.13  shows  the  growth  rate  vs. 
time  for  three  different  k  values  at  their  optimized  grid  setup.  These  curves  all  converge  to 
their  corresponding  constant  value. 

Figure  5.14  shows  the  comparison  of  computational  results  and  predicted  values  for 
linear  analysis.  It  shows  a  very  good  agreement  for  k  =  0.3, 0.5  and  0.7.  The  reason  for 
error  for  k  =  0.9  may  due  to  the  desire  of  even  fine  grid  for  the  computational  domain. 

The  linear  analysis  shows  that  the  predicted  growth  rate  for  a  liquid  column  without 
the  presence  of  gas  phase  are  matched  very  well  with  the  computational  results  from  3D 
BEM  model.  Thus,  the  validation  of  3D  BEM  model  is  completed. 
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Figure  5.2  Multiple  Values  of  Normal  Velocities  at  a  Corner  Node 
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Figure  5.4  Three  Different  Grid  Setups  for  a  Liquid  Jet 
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Figure  5.5  Schematic  Description  of  Surrounding  Grids  and  Elements  Around  a  Corner 

Node 
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Figure  5.7  Effect  of  Perturbation  Magnitude  on  Transverse  Modes  m  =  1  Mode:  k  =  1, 

Total  Length  =  3  A 


Figure  5.9  Effect  of  Perturbation  Magnitude  on  Transverse  Modes  m  =  3  Mode: 

Total  Length  =  3  A 


k 


Figure  5.10  Same  Grid  Set-Up  ( intJZ  =  10,  intJThe  =  10,  and  int-R  =  3)  for 
k  =  0.3, 0.7,  and  0.9.  Perturbation  Magnitude  is  0.01  for  m  =  0. 
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Figure  5.11  Track  of  Growth  Rate  with  Time  for  k  =  0.7,  Perturbation  Magnitude  =  0.01 
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Figure  5.13  Track  of  Growth  Rate  with  Time  for  Three  Different  k  values  at  Their 

Optimized  Grid  Setup 
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9  Appendix  E  -  Dense  Spray  Model  Description 


Heister,  S.  D.,  “Modeling  Primary  Atomization  Processes”,  AIAA  98-3837,  34th  AIAA 
Joint  Propulsion  Conference,  1998. 
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MODELING  PRIMARY  ATOMIZATION 

PROCESSES 

S.  D.  Heister  * 

January  28,  1999 


Abstract 

This  paper  summarizes  recent  progress  in  modeling 
complex  atomization  and  injector  internal  flow  pro¬ 
cesses  using  a  variety  of  numerical  techniques.  This 
program  has  been  developed  under  AFOSR  sponsor¬ 
ship  with  the  goal  of  characterizing  the  contributions 
of  atomization  processes  to  combustion  instabilities 
in  liquid  rocket  engines.  Complete  viscous,  unsteady 
simulations  of  spray  formation  (including  the  dense 
spray  region)  should  be  available  within  the  year. 

Introduction 

In  a  liquid  rocket  engine  (LRE)  combustion  cham¬ 
ber,  the  atomization  process  serves  as  a  precursor 
to  complex  vaporization,  mixing,  and  reaction  pro¬ 
cesses.  Not  only  is  the  atomization  process  important 
to  characterizing  the  steady-state  performance  of  a 
LRE,  but  it  also  can  play  an  important  role  in  un¬ 
steady  processes  within  the  combustion  chamber.  In 
fact,  numerous  authors1-5  have  implicated  the  atom¬ 
ization  process  as  a  mechanism  to  (at  least  partially) 
explain  high-frequency  LRE  combustion  instabilities. 

Even  under  steady  conditions,  our  knowledge  of 
the  detailed  processes  leading  to  atomization  of  a 
liquid  jet  is  quite  modest.  For  this  reason,  there 
have  been  relatively  few  efforts  aimed  at  improving 
our  understanding  of  injection  processes  under  dy¬ 
namic  conditions.  Most  previous  efforts  have  been 
experimental  in  nature  and  have  been  motivated  by 
LRE  combustion  stability  concerns.  Under  AFOSR 
sponsorship,  a  series  of  two-dimensional,  nonlinear 
models6-12  have  been  developed  to  ascertain  the  in¬ 
fluence  of  unsteady  chamber  conditions  on  basic  at¬ 
omization  processes.  These  models  utilize  Boundary 
Element  Methods  (BEMs)  to  provide  high  resolution 
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surface  shapes  approaching  and  beyond  actual  atom¬ 
ization  events. 

As  a  result  of  these  developments,  we  have  been 
able  to  quantify  the  dramatic  effects  chamber  un¬ 
steadiness  can  have  on  droplet  sizes  and  distortion 
of  liquid  jets/columns  emanating  from  impinging  ele¬ 
ment  injectors.  In  addition,  viscous  simulations  have 
revealed  that  a  “swelling”  of  the  free  surface  takes 
place  immediately  outside  the  orifice  exit.  This  phe¬ 
nomena,  which  has  been  observed  in  the  die-cast  in¬ 
dustry,  is  attributed  to  boundary  layer  rearrangement 
due  to  the  transition  from  a  no-slip  to  a  free-surface 
boundary  condition.  Results  indicate  the  degree  of 
swelling  is  directly  proportional  to  the  boundary  layer 
thickness  at  the  exit  of  the  orifice.  In  fact,  this  con¬ 
clusion  helps  us  understand  the  role  of  orifice  passage 
design  in  affecting  subsequent  atomization  processes. 

Recognizing  the  importance  of  the  internal  flow  to 
the  atomization  process,  we  have  also  developed  mod¬ 
els  capable  of  addressing  this  flowfield.  Efforts13-15 
in  this  area  have  focused  on  the  role  of  cavitation 
in  influencing  orifice  passage  flows.  At  the  present 
time,  we  have  a  fully  viscous,  unsteady,  3-D  capa¬ 
bility  to  resolve  either  single  or  two-phase  flows  in 
these  passageways.  While  cavitation  may  not  be  of 
critical  interest  in  high  pressure  rocket  engine  appli¬ 
cations,  the  two-phase  flow  methodology  is  useful  in 
extending  solutions  to  the  spray  region  formed  out¬ 
side  the  orifice.  The  approach  is  based  on  a  pseudo¬ 
density  treatment  in  which  the  pseudo-density  varies 
in  magnatude  between  the  liquid  and  vapor  (or  gas) 
extremes.  By  developing  a  constitutive  relation  for 
the  pseudo-density,  one  can  investigate  a  variety  of 
two- phase  flows. 

In  this  paper,  we  will  highlight  recent  developments 
of  a  3-D  BEM  tool  capable  of  solving  more  general 
atomization  problems.  In  addition,  we  will  briefly 
describe  the  pseudo-density  modeling  approach  and 
highlight  a  new  model  capable  of  resolving  the  entire 
spray  (including  dense  spray  region)  while  accounting 
for  gas  entrainment  and  momentum  transfer  between 
phases. 
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Boundary  Element  Modeling 

During  the  past  two  years,  we  have  been  developing 
a  fully  three-dimensional  capability  to  simulate  com¬ 
plex  fluid  surfaces  under  arbitrarily-large  distortions. 
Originally,  these  developments  were  undertaken  with 
Professor  Chandrajit  Bajaj  of  our  Computer  Science 
department.  Professor  Bajaj  Js  group  had  agreed  to 
implement  some  of  their  surface  tracking/regridding 
strategies,  while  our  group  focused  on  the  develop¬ 
ment  of  the  fluid  flow  solver.  Unfortunately,  Profes¬ 
sor  Bajaj  left  Purdue  for  a  position  at  the  University 
of  Texas  in  January  of  this  year;  for  this  reason,  we 
had  to  substantially  rescope  the  efforts. 

During  the  past  year,  we  have  developed  a  fully 
3-D,  unsteady  flow  solver  using  linear  Boundary  Ele¬ 
ments,  as  well  as  computational  mesh  generation  ca¬ 
pabilities  for  droplet  and  liquid  jet  geometries.  Fig¬ 
ure  1  highlights  computational  meshes  generated  for 
droplets  via  interpolation  from  octahedron  polygons. 
This  general  technique  can  also  be  applied  to  arbi¬ 
trary  initial  fluid/structural  shapes.  Using  the  oc¬ 
tahedron  grid  generation  package,  simulations  were 
performed  for  droplets  oscillating  in  2nd,  3rd,  and 
4th  modes.  An  example  simulation  result  is  shown 
in  Fig.  2  for  a  3rd  mode  oscillation.  Computed  fre¬ 
quencies  agree  well  with  those  in  the  literature  for 
linear/nonlinear  droplet  oscillations.  Further  valida¬ 
tion  was  conducted  on  the  droplet  case  to  insure  that 
solutions  were  independent  of  grid.  Furthermore,  a 
rather  exhaustive  study  was  completed  to  insure  con¬ 
vergence  of  all  the  complex  integrals  involved  in  the 
3-D  BEM. 

Having  validated  the  model,  we  set  out  to  study 
the  nonlinear  stability  of  an  infinite  liquid  jet.  The 
infinite-jet  assumption  leads  to  periodic  boundary 
conditions  on  a  computational  domain  involving  a 
single  wavelength  of  a  temporally-evolving  instabil¬ 
ity.  The  use  of  a  single  wavelength  reduces  the  size 
of  the  domain  dramatically,  but  introduces  substan¬ 
tial  complications  as  a  result  of  corners  created  at 
the  ends  of  the  wave.  A  typical  computational  mesh 
is  depicted  in  Fig.  3;  a  structured  triangular  mesh 
was  generated  in  this  case.  Using  this  mesh,  we  have 
evaluated  the  model  against  analytic  results  which 
predict  the  growth  rate  of  an  instability  of  a  given 
wavelength. 

Figure  4  provides  a  comparison  of  these  results 
for  axisymmetric  instabilities  imposed  on  an  axial 
wave  with  wavenumber  k .  The  solid  curve  represents 
the  analytic  result16,  while  the  BEM  calculations  are 
shown  as  symbols  on  the  plot.  Computational  re¬ 
sources  constrained  us  to  fairly  coarse  grids  including 
20-40  axial  nodes  and  10-20  circumferential  nodes  on 
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Figure  1:  Computational  Mesh  Generated  from  Oc¬ 
tahedron  Interpolations  for  use  in  Droplet  Oscillation 
Simulations 

the  interface.  As  a  result  of  these  coarse  meshes,  sub¬ 
stantial  disagreement  is  noted  between  model  results 
and  analytic  solutions.  In  addition,  the  results  point 
to  the  need  for  improved  surface  fitting  (our  current 
approach  provides  local  fitting  with  a  second-order 
polynomial  after  Chahine,  et.  al17).  Since  an  ad¬ 
vanced  surface  fitting  and  regridding  capability  has 
not  been  included  due  to  the  departure  of  Dr.  Ba¬ 
jaj,  the  present  model  does  not  have  the  capability 
to  handle  large  surface  distortions.  In  the  future,  we 
hope  to  develop  this  capability  ourselves. 

Another  area  to  be  explored  in  the  near  future  is 
the  simulation  of  a  pressure-swirl  atomizer  such  as 
those  used  in  numerous  cryogenic  and  Russian  rocket 
engines.  These  axisymmetric  simulations  will  be  ex¬ 
plored  through  modifications  of  existing  finite-length 
liquid  jet  codes.  Orifice  geometry,  swirl  velocity,  fluid 
properties,  and  general  unsteady  inflow  conditions 
can  effectively  be  considered  in  the  simulations.  We 
expect  to  have  simulations  of  this  flow  by  the  end  of 
1998. 

Viscous,  Two-Phase  Flow  Simu¬ 
lations 

There  is  a  critical  need  for  a  model  capable  of  address¬ 
ing  fully  coupled  spray  dynamics  in  which  momen- 
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Figure  2:  Droplet  Shapes  at  Various  Times;  3rd  Mode 
Oscillation 


Figure  3:  Computational  Mesh  for  Simulations  of 
Nonlinear  Stability  of  Infinite  Jet 
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Figure  4:  Replication  of  Analytic  Growth  Rates  Using 
the  Infinite  Jet  BEM .  Here ,  to  is  the  Growth  Rate ,  k  is 
the  Axial  Wavenumber ,  and  n  is  the  Circumferential 
Wavenumber 


turn  transfer  between  phases  is  handled  in  an  implicit 
fashion.  Models  based  on  linear  theory  are  known  to 
be  inaccurate,  yet  full  multidimensional  simulations 
are  well  beyond  computational  capabilities.  What 
we  really  desire  is  a  mechanism  to  obtain  averaged 
flowfield  conditions  which  are  a  reflection  of  the  lo¬ 
cal  droplet  number  density  (or  local  “void  fraction”) 
without  having  to  compute  the  flowfield  around  all 
the  individual  droplets.  Models  of  this  type  have 
been  in  use  for  many  years  in  assessing  two-phase 
bubbly  flows.  Recently,  models  of  this  type  have  been 
successfully  developed  to  assess  unsteady  behavior  in 
cavitating  flows.  The  models  are  based  on  a  pseudo¬ 
density  formulation,  and  are  sometimes  referred  to 
as  “single-fluid”  models.  Here,  the  pseudo-density  is 
a  fictitious  parameter  which  varies  between  the  liq¬ 
uid  and  gas  density  extremes,  i.e.  a  pseudo-density 
value  midway  between  liquid  and  gas  densities  would 
be  representative  of  a  mixture  which  is  50%  liquid 
and  50%  gaseous. 

Cavitating  Internal  Orifice  Flows 

The  single  fluid  models  are  powerful  in  that  they  are 
capable  of  performing  fully  viscous,  unsteady  calcu¬ 
lations  for  complex  two-phase  flows.  The  main  chal¬ 
lenge  in  using  this  methodology  stems  from  the  fact 
that  there  is  no  equation  of  state  which  defines  the 
fictitious  pseudo-density.  Therefore  a  constitutive  re¬ 
lation  for  this  parameter  must  be  developed  based  on 
the  characteristics  of  the  two-phase  flow.  In  cavitat¬ 
ing  flows,  we  have  been  successful  in  deriving  a  con- 
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Figure  5:  Schematic  of  2-D  Slot  Experimental 
Apparatus 18 

stitutive  relation  for  the  pseudo-density  based  on  the 
response  of  a  single  bubble  to  a  time- varying  external 
pressure  (Rayleigh-Plesset  equation).  Using  this  ap¬ 
proach,  we  have  developed  a  time-dependent  relation 
for  the  pseudo-density  which  is  reflective  of  local  flow- 
field  conditions  and  the  pressure  history  encountered 
by  a  given  bubble. 

Using  this  implementation,  we  have  successfully 
modeled  complex,  unsteady  cavitation  processes  in 
both  external  and  internal  flows.  Recent  simulations 
have  focused  on  a  2-D  slot  flow  in  order  to  validate 
computations  against  recent  experimental  data  ob¬ 
tained  by  Henry18.  Figure  5  provides  a  schematic 
representation  of  the  experimental  apparatus  used  to 
obtain  a  single,  unambiguous  cavitation  surface  off 
the  lip  of  a  single  slot. 

Figures  6  and  7  provide  some  sample  comparisons 
of  the  model  with  the  experimental  data.  In  Fig. 
6,  we  compare  the  overall  length  of  the  cavity  as 
a  function  of  the  imposed  pressure  drop  for  a  slot 


Figure  6:  Comparison  of  Cavitation  Length  with  Ex¬ 
perimental  Data  of  Henry 18 

with  L/D  =  10.7.  Since  the  flowfield  tends  to  be 
unsteady,  we  have  included  both  minimum  and  max¬ 
imum  lengths  obtained  in  the  numerical  simulation 
with  Henry’s  data  which  was  obtained  at  a  given  in¬ 
stant  in  time.  Overall,  the  agreement  is  quite  good 
considering  the  complex,  turbulent  cavitation  pro¬ 
cesses  which  are  evident  in  the  experimental  images. 

A  comparison  with  an  actual  image  is  shown  in 
Fig.  7.  Here,  pseudo-density  contours  of  0.24,  0.49, 
0.74  and  0.99  (from  inside  to  outside)  are  shown  in 
the  broad  white  lines  in  the  image.  Once  again,  the 
computed  flow  is  quasi-periodic,  so  the  comparison 
is  made  for  just  one  instant  in  this  process;  there  is 
no  way  to  correlate  the  timing  with  the  single  experi¬ 
mental  image.  In  spite  of  this  fact,  comparisons  using 
other  points  in  the  periodic  process  are  quite  similar. 
The  model  predicts  two  smaller  regions  of  cavitation 
while  the  experimental  data  reveals  a  single  large  im¬ 
age.  Overall,  the  extent  of  cavitation  is  predicted 
well.  Turbulence  effects,  which  are  not  modeled  in 
the  present  code,  are  a  probable  explanation  for  the 
observed  discrepancies. 

Future  Work  -  Dense  Spray  Model 

While  the  pseudo-density  formulation  can  be  em¬ 
ployed  when  the  gas  is  the  disperse  phase,  we  must 
provide  an  alternate  constitutive  relation  for  p  when 
the  liquid  is  the  disperse  phase  (as  in  the  case  of  a 
spray).  Since  p  is  simply  a  reflection  of  the  local  liq- 
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Figure  7:  Comparison  of  Pseudo-Density  Contours 
(in  white)  with  Experimental  Image  for  A  P=  38  psi, 
L/D  —  10.7.  Here ,  the  Flow  Direction  is  Upward  and 
Outermost  Density  Contour  Corresponds  to  p  —  0.99. 


uid  volume  in  a  given  sample  space,  we  could  provide 
a  straightforward  calculation  of  pseudo-density  if  we 
knew  the  location  and  size  of  all  the  droplets  in  the 
flow.  However,  since  there  are  too  many  droplets  to 
track  each  one  individually,  we  must  consider  a  subset 
of  the  total  droplet  field  with  which  to  perform  our 
calculation.  If  we  presume  that  the  spray  can  be  rep¬ 
resented  by  several  discrete  planes  of  droplets,  each 
with  identical  characteristics,  we  could  effectively 
compute  pseudo-density  by  considering  the  dynam¬ 
ics  of  the  droplets  lying  in  just  one  of  these  planes. 
This  procedure  amounts  to  the  presumption  that  the 
spray  is  “quasi- axisymmetric” ;  that  the  droplet  sizes 
and  velocities  are  similar  in  any  azimuthal  plane  in¬ 
tercepting  the  main  axis  of  the  spray.  Numerous  ex¬ 
periments  indicate  that  such  a  symmetry  exists  for 
nozzles,  especially  for  the  predominate  case  where  the 
orifice  passage  is  axisymmetric. 

The  unique  treatment  of  a  quasi- axisymmetric 
droplet  field  will  provide  the  capability  to  track  a 
manageable  number  of  droplets  in  a  single  plane. 
While  there  is  little  research  in  the  dense  sprays 
arena,  much  is  known  about  the  dynamics  of  sin¬ 
gle  droplets  and  groups  of  droplets.  In  recent 
years,  we  have  seen  advancements  in  knowledge  of 
droplet  collisions19’20,  secondary  atomization2 1~24, 
evaporation6'5*’  25~28  behavior  in  both  steady  and 
unsteady  environments,  and  a  wealth  of  knowledge 
regarding  group  combustion  and  dynamics.  It  makes 
good  sense  to  fold  this  knowledge  into  a  comprehen¬ 
sive  formulation  for  an  entire  spray.  The  recent  devel¬ 
opment  of  pseudo-density  codes  and  the  vast  research 
on  droplet  dynamics  provides  a  unique  opportunity 
for  the  development  of  the  first  comprehensive  spray 
model. 

Initial  Conditions 

Prescribing  a  set  of  well-posed  initial  conditions  is  al¬ 
ways  problematic  in  spray  modeling.  The  most  nat¬ 
ural  approach  would  simulate  the  opening  of  a  valve 
and  the  motion  of  the  initial  slug  of  fluid  through  the 
orifice  passage.  However,  this  option  is  not  viable 
using  the  current  methodology  since  the  complex  in¬ 
terface  tracking  would  introduce  substantial  compli¬ 
cations  to  the  envisioned  algorithm.  Moreover,  using 
this  approach,  one  would  have  to  perform  very  long 
(and  costly)  calculations  in  order  to  reach  a  quasi¬ 
steady  condition  in  the  spray. 

Given  this  state  of  affairs,  we  propose  the  geome¬ 
try  and  initial  conditions  highlighted  in  Fig.  3.  Ba¬ 
sic  inputs  to  the  model  include  fluid  properties,  ori¬ 
fice  diameter  (D),  injection  velocity,  (i>;),  jet  intact 
core  length,  (Lc),  spray  cone  angle,  (a),  and  initial 
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Figure  8:  Schematic  Denoting  Assumed  Ini¬ 

tial  Conditions  and  Coordinate  System  for  Quasi- 
Axisymmetric  Model.  Large  Droplets  of  Uniform 
Size  are  “ Launched ”  from  the  Conical  Surface  Shown . 
The  Cone  Angle,  a,  and  Intact  Core  Length,  Lc  are 
Input  from  Experimental  Observations 


droplet  diameter,  (c?2).  Droplets  are  “launched”  into 
the  computational  domain  along  the  conical  surface 
emanating  from  the  lip  of  the  orifice  and  terminating 
at  the  point  z  —  Lc,  r  —  0.  This  approach  is  required 
since  we  intend  to  search  for  droplet  collisions;  there 
is  insufficient  surface  area  available  at  the  exit  plane 
to  provide  the  proper  massflux  of  droplets  without 
having  these  droplets  overlap.  This  is  an  important 
distinction,  since  simulations  using  the  KIVA  code29 
launch  large  “blobs”  of  fluid  which  represent  large 
numbers  of  drops  and  overlap  to  a  great  degree  in 
order  to  simulate  droplet  injection  into  the  field. 

The  computational  domain  will  include  the  edge  of 
the  conical  surface  bounded  by  a  large  cylinder.  On 
the  outer  boundaries,  constant  pressure  conditions 
can  be  assumed  and  velocities  can  be  extrapolated 
from  the  interior  of  the  domain.  Individual  droplets 
in  the  spray  will  pass  through  the  “end”  of  the  do¬ 
main  at  a  location,  say,  2  ~  zmax.  Spray  statistics 
including  droplet  sizes  and  velocities  can  be  obtained 
at  this  boundary  for  direct  comparison  with  PDPA  or 
Malvern  experimental  data.  In  fact,  one  can  develop 
sets  of  spray  statistics  at  arbitrary  z  locations  using 
the  proposed  methodology. 

Since  the  initial  droplet  size,  cfc,  cannot  be  read¬ 
ily  determined  from  experimental  or  analytic  means, 
there  is  a  concern  that  this  input  may  drastically  ef¬ 
fect  the  results.  However,  if  we  launch  large  drops,  we 
can  assure  that  all  drops  go  through  numerous  sec¬ 
ondary  atomization  and  collision  processes  thereby 
making  results  less  sensitive  to  this  initial  input. 


Droplet  Dynamics  Module 

Probably  the  most  crucial  element  of  our  determin¬ 
istic  model  of  a  fully-coupled  dense  spray  is  the  re¬ 
sponse  of  individual  droplets  to  the  imposed  flowfield. 
At  the  present  time,  no  one  has  undertaken  the  ar¬ 
duous  task  of  describing  secondary  atomization  and 
collision  processes  for  a  dense  spray.  However,  we  be¬ 
lieve  the  time  is  right  for  such  a  development  for  the 
following  reasons: 

•  Recent  work  at  Princeton  provides  a  complete 
description  of  droplet  collision/coalescence  for 
both  on  and  off-axis  collisions  of  droplets  com¬ 
posed  of  various  fluids.  These  results  can  eas¬ 
ily  be  incorporated  parametrically  into  the  pro¬ 
posed  model. 

•  Recent  work  at  Michigan  and  Wisconsin  provides 
a  quantitative  description  of  secondary  atomiza¬ 
tion  process  in  terms  of  dimensionless  parame¬ 
ters  such  as  Weber  and  Ohnesorge  numbers.  Al¬ 
gebraic  equations  already  exist  for  drop  sizes  re¬ 
sulting  from  the  secondary  atomization  process, 
and  effective  drag  coefficient  data  is  also  avail¬ 
able  from  a  variety  of  sources. 

•  Advances  in  computational  power  and  parallel 
processing  architectures  make  it  possible  to  track 
tens,  if  not  hundreds  of  thousands  of  individual 
droplets  in  a  single  calculation.  While  this  alone 
is  still  not  adequate  for  simulation  of  a  complete 
spray,  the  unique  methodology  described  previ¬ 
ously  will  permit  a  quasi-axisymmetric  model  to 
be  developed. 

With  these  ideas  in  mind,  we  focus  on  a  single  drop 
at  some  arbitrary  location  within  the  spray  as  shown 
in  Fig.  4. 

We  presume  a  coordinate  system  moving  with  the 
drop  in  a  Lagrangian  fashion.  The  x  axis  is  aligned 
with  the  drop’s  velocity  vector  such  that  drag  forces 
act  along  this  axis  and  “lift”  forces  act  along  the  ac¬ 
companying  y  axis.  Note  that  a  lift  force  can  be 
present  since  the  pseudo-fluid  velocity  vector  will  gen¬ 
erally  not  be  aligned  with  the  droplet  velocity  vector 
(as  shown  in  Fig.  4).  This  situation  arises  due  to  the 
fact  that  a  group  of  droplets  is  present;  local  pseudo¬ 
density  gradients  manifested  by  the  presence  of  neigh¬ 
boring  droplets  will  account  for  droplet /wake  inter¬ 
actions  and  other  group  dynamic  effects.  Assuming 
the  drop  is  much  denser  than  the  surrounding  gas,  we 
can  effectively  neglect  unsteady  virtual  mass  and  Bas¬ 
set  forces  such  that  drag  and  lift  are  the  only  forces 
acting  on  the  drop.  Under  this  situation,  Newton’s 
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Figure  9:  Motion  of  a  Typical  Drop  within  the  Spray. 
Local  Properties  of  the  Pseudo-Fluid  (v )  are  Assumed 
to  be  Known  as  a  Result  of  Solution  of  Navier-Stokes 
Equations 


second  law  gives: 

Md-jjj^- =  CjyArtq  (1) 

Mdjt |  =  CLALq  (2) 

where  Md  is  the  mass  of  the  droplet,  Cd  and  Cl  are 
drag  and  lift  coefficients,  and  td  is  the  time  as  mea¬ 
sured  from  the  launching  point  for  a  given  droplet. 
Assuming  the  droplet  deforms  into  an  ellipsoid,  the 
effective  cross-sectional  areas  Ad  and  Al  can  be  es¬ 
tablished  using  the  initial  and  current  droplet  radii. 
Determination  of  the  effective  dynamic  pressure  act¬ 
ing  on  the  drop  ( q )  is  not  trivial  in  this  case  since  the 
velocity  returned  from  the  Navier  Stokes  solver  (iT)  is 
a  reflection  of  the  local  pseudo-density  and  not  the 
density  of  the  disperse  phase.  We  propose  to  initiate 
our  developments  assuming: 

[(N-l?l))IN-M)0  (3) 

This  form  effectively  treats  the  pseudo-fluid  veloc¬ 
ity,  v ,  as  that  of  a  disperse  phase,  but  it  does  provide 
the  correct  asymptotic  limits  for  both  dilute  and  very 
dense  regions  within  the  spray.  In  the  dilute  region 
of  the  spray  (p  — >  pg,  |tT|  -»  0),  and  the  drag  on 


the  drop  is  equivalent  to  a  single  drop  propagating 
into  a  quiescent  gas.  In  the  dense  region  of  the  spray 
(p  pu  v  ->  Vd  and  Eq.  6  reduces  to  q  =  0  so  that 
the  drag  vanishes  in  this  limit.  Note  that  the  as¬ 
sumed  form  for  q  incorporates  the  “blockage  effect” 
realized  when  a  droplet  is  traveling  in  the  wake  of  one 
or  more  ofter  droplets.  Since  there  is  little  research 
available  to  draw  upon  for  further  insight,  the  rele¬ 
vancy  of  Eq.  6  will  have  to  be  demonstrated  as  a  part 
of  the  model  development.  Using  this  methodology 
and  empirical  data  for  force  coefficients,  one  can  in¬ 
tegrate  the  trajectory  of  any  droplet  in  the  flow.  To 
track  the  droplet  deformation  with  time,  we  propose 
to  use  the  technique  of  Ibrahim,  et.  al49.  Knowl¬ 
edge  of  the  droplet  deformation  can  be  used  to  set 
breakup  criteria  in  the  secondary  atomization  model 
described  below. 


Droplet  Collision/Secondary  Atomization 

Droplet /droplet  collisions  are  known  to  be  an  impor¬ 
tant  physical  process  occurring  within  a  spray.  Since 
droplets  in  the  dilute  region  are  decelerated  due  to  the 
higher  drag  with  the  gaseous  phase,  they  are  continu¬ 
ally  subject  to  collisions  with  faster  moving  droplets 
emanating  from  the  interior  of  the  spray.  Impor¬ 
tant  momentum  exchange  and  secondary  atomization 
takes  place  as  a  result  of  this  process. 

Due  to  the  recent  basic  research  efforts  in  Professor 
Law’s  group  at  Princeton37,38,  we  are  in  a  good  posi¬ 
tion  to  account  for  these  physical  processes  in  the  pro¬ 
posed  model.  The  Princeton  research  has  shown  that 
the  collision  process  can  be  characterized  in  terms 
of  a  collision  parameter,  J5,  and  the  effective  Weber 
number  based  on  the  relative  collision  velocity,  Wec: 

B  =  2 X/{d1+d2)  Wec  =  p,(d1  +  d2)v2c/(2cr)  (4) 

where  X  is  distance  between  droplet  centers  at  the 
impact  point,  d\  and  d%  are  the  diameters  of  the  col¬ 
liding  drops,  and  vc  is  the  relative  collision  velocity. 
The  type  of  liquid  and  overall  pressure  level  have  also 
been  shown  to  be  important  in  the  collision  process; 
at  higher  pressures  the  results  appear  to  be  relatively 
insensitive  to  this  parameter. 

Both  virgin  droplets  and  those  formed  as  a  re¬ 
sult  of  a  collision  process  are  subject  to  secondary 
atomization  caused  by  interaction  with  the  gaseous 
phase.  Within  the  past  few  years,  much  quantitative 
information  regarding  this  process  has  been  devel¬ 
oped,  primarily  by  Professor  Faeth’s  research  group 
at  Michigan.  Results  from  these  studies39,40  provide 
predicted  droplet  sizes  as  a  function  of  the  size  of  the 
parent  drop  and  the  effective  Weber  and  Ohnesorge 
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numbers  characterizing  the  initial  aerodynamic  con¬ 
ditions  to  which  the  droplet  was  subjected. 

There  are  two  problems  in  using  these  data  directly. 
Since  the  experiments  were  conducted  by  exposing 
the  droplet  to  a  uniform  flow  behind  a  weak  shock 
wave,  the  droplet  senses  a  dynamic  pressure  which 
is  a  maximum  at  the  initial  stages  of  the  trajectory. 
In  a  dense  spray,  the  droplet  is  shielded  such  that 
the  dynamic  pressure  may  be  an  increasing  function 
of  time.  Therefore,  it  is  not  immediately  clear  how 
to  characterize  the  Weber  and  Ohnesorge  numbers 
under  these  conditions.  This  issue  will  need  to  be 
resolved  through  a  series  of  numerical  experiments 
in  which  effective  Weber  and  Ohnesorge  numbers  are 
utilized  to  characterize  the  process. 

The  second  problem  arises  from  the  breakup  cri¬ 
teria  utilized  in  the  experiments.  Many  previous  ex¬ 
periments  using  single  droplets  exposed  to  decaying 
dynamic  pressure  histories  show  that  breakup  occurs 
at  a  specific  dimensionless  time.  Since  the  droplets 
emanating  from  the  spray  are  subjected  to  a  differ¬ 
ent  dynamic  pressure  history,  we  propose  to  utilize  a 
breakup  criteria  based  on  the  overall  deformation  of 
the  drop. 


Pseudo-Density  Routine 

Through  the  use  of  droplet  dynamics  module,  we 
should  be  able  to  determine  the  droplet  size  distribu¬ 
tion  throughout  space  (r,  z)  at  any  instant  in  time. 
By  determining  the  liquid  volume  in  any  fixed  sample 
volume,  Vs,  we  can  directly  calculate  p : 


P=Vi/Vs  (5) 

where  VJ  is  the  total  liquid  volume  lying  within  the 
selected  sample  volume.  The  main  challenge  lies  in 
picking  a  sample  volume  which  is  small  enough  to 
reflect  local  flowfield  conditions,  yet  large  enough  to 
permit  p  to  vary  smoothly. 

Figure  6  describes  the  methodology  for  determin¬ 
ing  the  local  pseudo-density  given  at  sample  radius, 
Rs.  The  quantity  Rdi  is  the  distance  to  droplet  “i” 
as  measured  from  the  center  of  the  sample  volume. 
We  must  take  into  account  the  three-dimensionality 
introduced  by  the  droplet  field  since  the  droplet  vol¬ 
ume  lies  outside  the  infinitely-thin  computational 
plane.  Moreover,  we  must  consider  droplets  in  ad¬ 
jacent  planes  which  enter  the  assumed  sample  vol¬ 
ume  in  order  to  get  an  accurate  reflection  of  the  local 
pseudo-density. 

If  we  choose  an  Rs  value,  then  we  can  evaluate  p  by 
adding  up  all  droplet  volumes  which  lie  in  the  sphere 


Figure  10:  Determining  the  Local  Pseudo-Density  by 
Searching  for  Droplets  which  are  Included  in  the  As¬ 
sumed  Sample  Volume.  Droplets  in  Adjacent  Planes 
also  must  be  Considered  in  this  Calculation. 


of  influence: 


Rdi  <  Rs  (6) 


In  general,  there  will  be  droplets  which  intersect  the 
boundary  Rs ,  but  we  will  neglect  this  factor  by  ac¬ 
cepting  only  droplets  whose  center  lies  within  the 
stated  boundary.  On  average  this  approach  is  jus¬ 
tified  since  it  will  be  equally  likely  for  the  droplet 
center  to  lie  on  either  side  of  this  boundary. 

At  the  present  time,  we  have  completed  prelimi¬ 
nary  versions  of  all  the  elements  of  the  droplet  dy¬ 
namics  module  and  have  a  capability  to  track  of  the 
order  of  105  droplets  in  the  computational  plane.  In 
addition,  droplet  collision/coalescence  and  secondary 
atomization  submodels  are  nearly  complete.  A  sam¬ 
ple  calculation  for  droplets  entering  a  vaccuum  (i.e. 
no  drag  model)  is  shown  in  Fig.  11.  Here,  the  evo¬ 
lution  of  the  spray  is  detailed  in  a  series  of  images 
at  different  times  from  the  injection  event.  When 
complete,  we  will  be  able  to  build  up  databases  of 
complete  spray  statistics  (drop  size  and  velocity  dis¬ 
tributions)  at  arbitrary  axial  or  radial  stations. 

We  anticipate  that  the  model  will  be  integrated 
with  the  viscous  flow  solver  sometime  this  fall  and 
checkout  of  test  cases  will  begin  in  early  1999.  When 
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Figure  11:  Sample  Result  of  Dense  Spray  Model  De¬ 
picting  Droplet  Field  Developing  in  Time . 


complete,  the  model  will  be  capable  of  addressing  a 
wide  array  of  spray  problems  due  to  its  general  na¬ 
ture. 


Conclusions 

This  paper  has  summarized  both  current  and  future 
efforts  in  the  modeling  of  primary  atomization  pro¬ 
cesses  with  emphasis  on  liquid  rocket  engine  applica¬ 
tions.  The  3-D  boundary  element  model  is  complete, 
but  will  require  surface  fitting  and  regridding  mod¬ 
ules  to  become  a  truly  effective  tool.  Models  based 
on  a  pseudo-density  formulation  of  the  Navier- Stokes 
equations  show  much  promise  for  these  complex  two- 
phase  flows.  We  have  successfully  used  such  a  model 
to  describe  viscous,  unsteady,  cavitating  internal  in¬ 
jector  flows.  Using  similar  notions,  we  are  currently 
developing  the  first  model  capable  of  investigating  the 
fully  coupled  dynamics  of  a  dense  spray  under  both 
quasi-steady  and  unsteady  conditions.  This  general 
tool  should  be  of  great  interest  to  the  atomization 
community  when  completed  and  validated  against  ex¬ 
perimental  results. 
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